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

Last change on this file since 31466 was 31466, checked in by sjamaan, 6 years ago

numbers: Convert quotient to new style

File size: 81.1 KB
Line 
1/* numbers-c.c
2 *
3 * Copyright 2010-2014 The CHICKEN Team
4 *
5 * This contains a barely recognizable version of c/bignum.c from Scheme48 1.8:
6 * Copyright (c) 1993-2008 Richard Kelsey and Jonathan Rees
7 * Copyright 1986,1987,1988,1989,1990,1991 Massachusetts Institute of Technology
8 * Copyright 1992,1993,1994,2004 Massachusetts Institute of Technology
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions are
12 * met:
13 *
14 *    1. Redistributions of source code must retain the above copyright
15 *       notice, this list of conditions and the following disclaimer.
16 *
17 *    2. Redistributions in binary form must reproduce the above
18 *       copyright notice, this list of conditions and the following
19 *       disclaimer in the documentation and/or other materials provided
20 *       with the distribution.
21 *
22 *    3. The name of the author may not be used to endorse or promote
23 *       products derived from this software without specific prior
24 *       written permission.
25 *
26 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
27 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
28 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
29 * DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
30 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
31 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
32 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
33 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
34 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
35 * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
36 * POSSIBILITY OF SUCH DAMAGE.
37 */
38
39#include <assert.h>
40#include <errno.h>
41#include <math.h> /* frexp() */
42
43#define nmax(x, y)     ((x) > (y) ? (x) : (y)) /* From runtime.c */
44
45static void *tags;
46
47#include "numbers-c.h"
48
49static C_word init_tags(___scheme_value tagvec);
50static C_word check_number(C_word x);
51static void bignum_negate_2(C_word c, C_word self, C_word new_big);
52static void allocate_bignum_2(C_word c, C_word self, C_word bigvec);
53static C_word bignum_digits_destructive_scale_up_with_carry(C_word *start, C_word *end, C_word fix_factor, C_word carry);
54
55static void bignum_plus_unsigned(C_word k, C_word x, C_word y, C_word negp);
56static void bignum_plus_unsigned_2(C_word c, C_word self, C_word result);
57static int bignum_cmp_unsigned(C_word x, C_word y);
58static void bignum_minus_unsigned(C_word k, C_word x, C_word y);
59static void bignum_minus_unsigned_2(C_word c, C_word self, C_word result);
60
61static void integer_times_2(C_word c, C_word self, C_word new_big);
62static C_regparm void bignum_digits_multiply(C_word x, C_word y, C_word result);
63static void bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp);
64static void bignum_times_bignum_unsigned_2(C_word c, C_word self, C_word result);
65static void digits_to_integer_2(C_word c, C_word self, C_word result);
66static void bignum_to_digits_2(C_word c, C_word self, C_word working_copy);
67static void bignum_to_digits_3(C_word c, C_word self, C_word string);
68static void flo_to_int_2(C_word c, C_word self, C_word result);
69static C_word ilen(C_uword x);
70static void bignum_allocate_for_shift(C_word c, C_word self, C_word x);
71static void bignum_negate_after_shift(C_word c, C_word self, C_word result);
72static void bignum_actual_shift(C_word c, C_word self, C_word result);
73static void bignum_abs_2(C_word c, C_word self, C_word result);
74static void bignum_random_2(C_word c, C_word self, C_word result);
75static void bignum_maybe_negate_magnitude(C_word k, C_word result);
76static void bignum_divrem_fixnum_2(C_word c, C_word self, C_word negated_big);
77static void bignum_negneg_bitwise_op(C_word c, C_word self, C_word result);
78static void bignum_posneg_bitwise_op(C_word c, C_word self, C_word result);
79static void bignum_pospos_bitwise_op(C_word c, C_word self, C_word result);
80static void C_ccall quotient_intflo_2(C_word c, C_word self, C_word intnum);
81static void C_ccall quotient_intflo_3(C_word c, C_word self, C_word intnum);
82static void bignum_destructive_normalize(C_word target, C_word source, C_word shift_left);
83static void bignum_destructive_remainder_unsigned_halfdigit(C_word k, C_word n, C_word d, C_word negp);
84static void bignum_destructive_divide_unsigned_halfdigit(C_word c, C_word self, C_word quotient);
85static void bignum_destructive_divide_unsigned_digit(C_word c, C_word self, C_word quotient);
86static C_word bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q);
87static C_word bignum_divide_and_subtract_digit(C_word v1, C_word v2, C_word guess, C_word *u);
88static void bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient);
89static void bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder);
90static void bignum_divide_2_unsigned_3(C_word c, C_word self, C_word tmp_big);
91static void bignum_destructive_divide_normalized(C_word u, C_word v, C_word q);
92static C_word bignum_divide_and_subtract(C_word *v_start, C_word *v_end, C_word guess, C_word *u_start);
93static C_word bignum_normalize_shifted(C_word bignum, C_word shift_right);
94
95static C_word
96init_tags(___scheme_value tagvec)
97{
98  tags = CHICKEN_new_gc_root();
99  CHICKEN_gc_root_set(tags, tagvec);
100  return C_SCHEME_UNDEFINED;
101}
102
103static C_word
104check_number(C_word x)
105{
106  C_word tagvec;
107
108  if((x & C_FIXNUM_BIT) != 0) return C_fix(FIX);
109  else if(C_immediatep(x)) return C_fix(NONE);
110  else if(C_header_bits(x) == C_FLONUM_TYPE) return C_fix(FLO);
111  else if(C_header_bits(x) == C_STRUCTURE_TYPE) {
112    tagvec = CHICKEN_gc_root_ref(tags);
113   
114    if (C_block_item(x, 0) == C_block_item(tagvec, BIG_TAG))
115      return C_fix(BIG);
116    else if (C_block_item(x, 0) == C_block_item(tagvec, COMP_TAG))
117      return C_fix(COMP);
118    else if (C_block_item(x, 0) == C_block_item(tagvec, RAT_TAG))
119      return C_fix(RAT);
120    else
121      return C_fix(NONE);
122  } else
123      return C_fix(NONE);
124}
125
126/* Copy all the digits from source to target, obliterating what was
127 * there.  If target is larger than source, the most significant
128 * digits will remain untouched.
129 */
130C_inline void bignum_digits_destructive_copy(C_word target, C_word source)
131{
132  C_memcpy(C_bignum_digits(target), C_bignum_digits(source),
133           /* TODO: This is currently in bytes.  If we change the
134            * representation that needs to change!
135            * We subtract the size of the header, too.
136            */
137           C_header_size(C_internal_bignum(source))-C_wordstobytes(1));
138}
139
140void C_ccall
141C_2_basic_plus(C_word c, C_word self, C_word k, C_word x, C_word y)
142{
143  C_word ab[nmax(C_SIZEOF_FIX_BIGNUM, C_SIZEOF_FLONUM*2)], *a = ab;
144
145  if (x & C_FIXNUM_BIT) {
146    if (y & C_FIXNUM_BIT) {
147      C_kontinue(k, C_a_u_i_2_fixnum_plus(&a, 2, x, y));
148    } else if (C_immediatep(y)) {
149      C_kontinue(k, C_SCHEME_FALSE);
150    } else if (C_block_header(y) == C_FLONUM_TAG) {
151      C_kontinue(k, C_flonum(&a, (double)C_unfix(x) + C_flonum_magnitude(y)));
152    } else if (C_IS_BIGNUM_TYPE(y)) {
153      C_u_2_bignum_plus(2, (C_word)NULL, k, C_a_u_i_fix_to_big(&a, x), y);
154    } else {
155      C_kontinue(k, C_SCHEME_FALSE);
156    }
157  } else if (C_immediatep(x)) {
158    C_kontinue(k, C_SCHEME_FALSE);
159  } else if (C_block_header(x) == C_FLONUM_TAG) {
160    if (y & C_FIXNUM_BIT) {
161      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) + (double)C_unfix(y)));
162    } else if (C_immediatep(y)) {
163      C_kontinue(k, C_SCHEME_FALSE);
164    } else if (C_block_header(y) == C_FLONUM_TAG) {
165      C_kontinue(k, C_a_i_flonum_plus(&a, 2, x, y));
166    } else if (C_IS_BIGNUM_TYPE(y)) {
167      C_kontinue(k, C_a_i_flonum_plus(&a, 2, x, C_a_u_i_big_to_flo(&a, 1, y)));
168    } else {
169      C_kontinue(k, C_SCHEME_FALSE);
170    }
171  } else if (C_IS_BIGNUM_TYPE(x)) {
172    if (y & C_FIXNUM_BIT) {
173      C_u_2_bignum_plus(2, (C_word)NULL, k, x, C_a_u_i_fix_to_big(&a, y));
174    } else if (C_immediatep(y)) {
175      C_kontinue(k, C_SCHEME_FALSE);
176    } else if (C_block_header(y) == C_FLONUM_TAG) {
177      C_kontinue(k, C_a_i_flonum_plus(&a, 2, C_a_u_i_big_to_flo(&a, 1, x), y));
178    } else if (C_IS_BIGNUM_TYPE(y)) {
179      C_u_2_bignum_plus(2, (C_word)NULL, k, x, y);
180    } else {
181      C_kontinue(k, C_SCHEME_FALSE);
182    }
183  } else {
184    C_kontinue(k, C_SCHEME_FALSE);
185  }
186}
187
188void C_ccall
189C_u_2_integer_plus(C_word c, C_word self, C_word k, C_word x, C_word y)
190{
191  C_word ab[C_SIZEOF_FIX_BIGNUM * 2], *a = ab;
192
193  if (x & C_FIXNUM_BIT && y & C_FIXNUM_BIT)
194    C_kontinue(k, C_a_u_i_2_fixnum_plus(&a, 2, x, y));
195
196  if (x & C_FIXNUM_BIT)
197    x = C_a_u_i_fix_to_big(&a, x);
198  if (y & C_FIXNUM_BIT)
199    y = C_a_u_i_fix_to_big(&a, y);
200
201  C_u_2_bignum_plus(2, (C_word)NULL, k, x, y);
202}
203
204/* Needs C_SIZEOF_FIX_BIGNUM */
205C_regparm C_word C_fcall
206C_a_u_i_2_fixnum_plus(C_word **ptr, C_word n, C_word x, C_word y)
207{
208  /* Exceptional situation: this will cause a real overflow */
209  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) && y == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
210    return C_bignum2(ptr, 1, 0, 2);
211  } else {
212    C_word z = C_unfix(x) + C_unfix(y);
213
214    if(!C_fitsinfixnump(z)) {
215      /* TODO: function/macro returning either fixnum or bignum from a C int */
216      /* This should help with the C API/FFI too. */
217      return C_bignum2(ptr, (z < 0), labs(z) & (C_uword)C_BIGNUM_DIGIT_MASK, 1);
218    } else {
219      return C_fix(z);
220    }
221  }
222}
223
224void C_ccall
225C_u_2_bignum_plus(C_word c, C_word self, C_word k, C_word x, C_word y)
226{
227  if (C_bignum_negativep(x)) {
228    if (C_bignum_negativep(y)) {
229      bignum_plus_unsigned(k, x, y, C_SCHEME_TRUE);
230    } else {
231      bignum_minus_unsigned(k, y, x);
232    }
233  } else {
234    if (C_bignum_negativep(y)) {
235      bignum_minus_unsigned(k, x, y);
236    } else {
237      bignum_plus_unsigned(k, x, y, C_SCHEME_FALSE);
238    }
239  }
240}
241
242void C_ccall
243C_2_basic_minus(C_word c, C_word self, C_word k, C_word x, C_word y)
244{
245  C_word ab[nmax(C_SIZEOF_FIX_BIGNUM, C_SIZEOF_FLONUM*2)], *a = ab;
246
247  if (x & C_FIXNUM_BIT) {
248    if (y & C_FIXNUM_BIT) {
249      C_kontinue(k, C_a_u_i_2_fixnum_minus(&a, 2, x, y));
250    } else if (C_immediatep(y)) {
251      C_kontinue(k, C_SCHEME_FALSE);
252    } else if (C_block_header(y) == C_FLONUM_TAG) {
253      C_kontinue(k, C_flonum(&a, (double)C_unfix(x) - C_flonum_magnitude(y)));
254    } else if (C_IS_BIGNUM_TYPE(y)) {
255      C_u_2_bignum_minus(2, (C_word)NULL, k, C_a_u_i_fix_to_big(&a, x), y);
256    } else {
257      C_kontinue(k, C_SCHEME_FALSE);
258    }
259  } else if (C_immediatep(x)) {
260    C_kontinue(k, C_SCHEME_FALSE);
261  } else if (C_block_header(x) == C_FLONUM_TAG) {
262    if (y & C_FIXNUM_BIT) {
263      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) - (double)C_unfix(y)));
264    } else if (C_immediatep(y)) {
265      C_kontinue(k, C_SCHEME_FALSE);
266    } else if (C_block_header(y) == C_FLONUM_TAG) {
267      C_kontinue(k, C_a_i_flonum_difference(&a, 2, x, y)); /* XXX NAMING! */
268    } else if (C_IS_BIGNUM_TYPE(y)) {
269      C_kontinue(k, C_a_i_flonum_difference(&a, 2, x, C_a_u_i_big_to_flo(&a, 1, y)));
270    } else {
271      C_kontinue(k, C_SCHEME_FALSE);
272    }
273  } else if (C_IS_BIGNUM_TYPE(x)) {
274    if (y & C_FIXNUM_BIT) {
275      C_u_2_bignum_minus(2, (C_word)NULL, k, x, C_a_u_i_fix_to_big(&a, y));
276    } else if (C_immediatep(y)) {
277      C_kontinue(k, C_SCHEME_FALSE);
278    } else if (C_block_header(y) == C_FLONUM_TAG) {
279      C_kontinue(k, C_a_i_flonum_difference(&a, 2, C_a_u_i_big_to_flo(&a, 1, x), y));
280    } else if (C_IS_BIGNUM_TYPE(y)) {
281      C_u_2_bignum_minus(2, (C_word)NULL, k, x, y);
282    } else {
283      C_kontinue(k, C_SCHEME_FALSE);
284    }
285  } else {
286    C_kontinue(k, C_SCHEME_FALSE);
287  }
288}
289
290void C_ccall
291C_u_2_integer_minus(C_word c, C_word self, C_word k, C_word x, C_word y)
292{
293  C_word ab[C_SIZEOF_FIX_BIGNUM * 2], *a = ab;
294
295  if (x & C_FIXNUM_BIT && y & C_FIXNUM_BIT)
296    C_kontinue(k, C_a_u_i_2_fixnum_minus(&a, 2, x, y));
297
298  if (x & C_FIXNUM_BIT)
299    x = C_a_u_i_fix_to_big(&a, x);
300  if (y & C_FIXNUM_BIT)
301    y = C_a_u_i_fix_to_big(&a, y);
302
303  C_u_2_bignum_minus(2, (C_word)NULL, k, x, y);
304}
305
306/* Needs C_SIZEOF_FIX_BIGNUM */
307C_regparm C_word C_fcall
308C_a_u_i_2_fixnum_minus(C_word **ptr, C_word n, C_word x, C_word y)
309{
310  C_word z = C_unfix(x) - C_unfix(y);
311
312  if(!C_fitsinfixnump(z)) {
313    /* TODO: function/macro returning either fixnum or bignum from a C int */
314    /* This should help with the C API/FFI too. */
315    return C_bignum2(ptr, (z < 0), labs(z) & (C_uword)C_BIGNUM_DIGIT_MASK, 1);
316  } else {
317    return C_fix(z);
318  }
319}
320
321void C_ccall
322C_u_2_bignum_minus(C_word c, C_word self, C_word k, C_word x, C_word y)
323{
324  if (C_bignum_negativep(x)) {
325    if (C_bignum_negativep(y)) {
326      bignum_minus_unsigned(k, y, x);
327    } else {
328      bignum_plus_unsigned(k, x, y, C_SCHEME_TRUE);
329    }
330  } else {
331    if (C_bignum_negativep(y)) {
332      bignum_plus_unsigned(k, x, y, C_SCHEME_FALSE);
333    } else {
334      bignum_minus_unsigned(k, x, y);
335    }
336  }
337}
338
339static void
340bignum_plus_unsigned(C_word k, C_word x, C_word y, C_word negp)
341{
342  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
343
344  if (C_bignum_size(y) > C_bignum_size(x)) {
345    C_word z = x;
346    x = y;
347    y = z;
348  }
349
350  k2 = C_closure(&ka, 4, (C_word)bignum_plus_unsigned_2, k, x, y);
351 
352  size = C_fix(C_bignum_size(x) + 1); /* One more digit, for possible carry. */
353  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
354}
355
356static void
357bignum_plus_unsigned_2(C_word c, C_word self, C_word result)
358{
359  C_word k = C_block_item(self, 1),
360         x = C_block_item(self, 2),
361         y = C_block_item(self, 3),
362         *scan_x = C_bignum_digits(x),
363         *end_x = scan_x + C_bignum_size(x),
364         *scan_y = C_bignum_digits(y),
365         *end_y = scan_y + C_bignum_size(y),
366         *scan_r = C_bignum_digits(result),
367         *end_r = scan_r + C_bignum_size(result),
368         sum, carry = 0;
369
370  /* Move over the two numbers simultaneously, adding digits w/ carry. */
371  while (scan_y < end_y) {
372    sum = (*scan_x++) + (*scan_y++) + carry;
373    if (C_fitsinbignumdigitp(sum)) {
374      (*scan_r++) = sum;
375      carry = 0;
376    } else {
377      (*scan_r++) = sum & C_BIGNUM_DIGIT_MASK;
378      carry = 1;
379    }
380  }
381 
382  /* The end of y, the smaller number.  Propagate carry into the rest of x. */
383  if (carry) {
384    while (scan_x < end_x) {
385      sum = (*scan_x++) + 1;
386      if (C_fitsinbignumdigitp(sum)) {
387        (*scan_r++) = sum;
388        carry = 0;
389        break;
390      } else {
391        (*scan_r++) = sum & C_BIGNUM_DIGIT_MASK;
392      }
393    }
394  }
395
396  if (carry) { /* We must've reached the end of x, with still a carry. */
397    (*scan_r) = 1; /* No need to trim: we're now using the extra digit. */
398    /* No need to normalize: We started with a bignum and it only grew. */
399    C_kontinue(k, result);
400  } else {
401    /* No more carry.  Copy remaining part of x (if any) to result. */
402    while (scan_x < end_x)
403      (*scan_r++) = (*scan_x++);
404
405    assert(scan_r == end_r - 1);
406
407    *scan_r = 0; /* Ensure trimming works.  TODO: set length directly? */
408    C_kontinue(k, C_bignum_normalize(result));
409  }
410}
411
412static int
413bignum_cmp_unsigned(C_word x, C_word y)
414{
415  C_word xlen = C_bignum_size(x), ylen = C_bignum_size(y);
416
417  if (xlen < ylen) {
418    return -1;
419  } else if (xlen > ylen) {
420    return 1;
421  } else if (x == y) {
422    return 0;
423  } else {
424    C_word *startx = C_bignum_digits(x);
425    C_word *scanx = startx + xlen;
426    C_word *scany = C_bignum_digits(y) + ylen;
427
428    while (startx < scanx) {
429      C_word xdigit = (*--scanx);
430      C_word ydigit = (*--scany);
431      if (xdigit < ydigit)
432        return -1;
433      if (xdigit > ydigit)
434        return 1;
435    }
436    return 0;
437  }
438}
439
440static void
441bignum_minus_unsigned(C_word k, C_word x, C_word y)
442{
443  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size, negp;
444
445  switch(bignum_cmp_unsigned(x, y)) {
446  case 0:             /* x = y, return 0 */
447    C_kontinue(k, C_fix(0));
448  case -1:            /* abs(x) < abs(y), return -(abs(y) - abs(x)) */
449    k2 = C_closure(&ka, 4, (C_word)bignum_minus_unsigned_2, k, y, x);
450   
451    size = C_fix(C_bignum_size(y)); /* Maximum size of result is length of y. */
452    C_allocate_bignum(3, (C_word)NULL, k2, size, C_SCHEME_TRUE, C_SCHEME_FALSE);
453  case 1:             /* abs(x) > abs(y), return abs(x) - abs(y) */
454    k2 = C_closure(&ka, 4, (C_word)bignum_minus_unsigned_2, k, x, y);
455   
456    size = C_fix(C_bignum_size(x)); /* Maximum size of result is length of x. */
457    C_allocate_bignum(3, (C_word)NULL, k2, size, C_SCHEME_FALSE, C_SCHEME_FALSE);
458    break;
459  }
460}
461
462static void
463bignum_minus_unsigned_2(C_word c, C_word self, C_word result)
464{
465  C_word k = C_block_item(self, 1),
466         x = C_block_item(self, 2),
467         y = C_block_item(self, 3),
468         *scan_x = C_bignum_digits(x),
469         *end_x = scan_x + C_bignum_size(x),
470         *scan_r = C_bignum_digits(result),
471         *scan_y = C_bignum_digits(y),
472         *end_y = scan_y + C_bignum_size(y),
473         difference, borrow = 0;
474
475  /* Move over the two numbers simultaneously, subtracting digits w/ borrow. */
476  while (scan_y < end_y) {
477    difference = (((*scan_x++) - (*scan_y++)) - borrow);
478    if (difference < 0) {
479      (*scan_r++) = ((C_word)1 << C_BIGNUM_DIGIT_LENGTH) + difference;
480      borrow = 1;
481    } else {
482      (*scan_r++) = difference;
483      borrow = 0;
484    }
485  }
486
487  /* The end of y, the smaller number.  Propagate borrow into the rest of x. */
488  if (borrow != 0) {
489    while (scan_x < end_x) {
490      difference = ((*scan_x++) - borrow);
491      if (difference < 0) {
492        (*scan_r++) = ((C_word)1 << C_BIGNUM_DIGIT_LENGTH) + difference;
493      } else {
494        (*scan_r++) = difference;
495        borrow = 0;
496        break;
497      }
498    }
499  }
500
501  assert(borrow == 0);
502
503  /* Finishing up: Copy remaining part of x (if any) into result. */
504  while (scan_x < end_x)
505    (*scan_r++) = (*scan_x++);
506
507  C_kontinue(k, C_bignum_normalize(result));
508}
509
510C_word C_ccall
511C_u_i_2_fixnum_gcd(C_word x, C_word y)
512{
513   C_word r;
514   
515   x = C_unfix(x);
516   y = C_unfix(y);
517   
518   if (x < 0) x = -x;
519   if (y < 0) y = -y;
520   
521   while(y != 0) {
522     r = x % y;
523     x = y;
524     y = r;
525   }
526   return C_fix(x);
527}
528
529C_word C_ccall
530C_a_u_i_2_flonum_gcd(C_word **p, C_word n, C_word x, C_word y)
531{
532   double xub, yub, r;
533   
534   xub = C_flonum_magnitude(x);
535   yub = C_flonum_magnitude(y);
536
537   if (xub < 0.0) xub = -xub;
538   if (yub < 0.0) yub = -yub;
539   
540   while(yub != 0.0) {
541     r = fmod(xub, yub);
542     xub = yub;
543     yub = r;
544   }
545   C_return(C_flonum(p, xub));
546}
547
548void C_ccall
549C_basic_abs(C_word c, C_word self, C_word k, C_word x)
550{
551  C_word ab[nmax(C_SIZEOF_FIX_BIGNUM, C_SIZEOF_FLONUM)], *a = ab;
552
553  if (x & C_FIXNUM_BIT)
554    C_kontinue(k, C_a_u_i_fixnum_abs(&a, 1, x));
555  else if (C_immediatep(x))
556    C_kontinue(k, C_SCHEME_FALSE);
557  else if (C_block_header(x) == C_FLONUM_TAG)
558    C_kontinue(k, C_a_i_flonum_abs(&a, 1, x));
559  else if (C_IS_BIGNUM_TYPE(x))
560    C_u_bignum_abs(1, (C_word)NULL, k, x);
561  else
562    C_kontinue(k, C_SCHEME_FALSE);
563}
564
565void C_ccall
566C_u_integer_abs(C_word c, C_word self, C_word k, C_word x)
567{
568  C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
569
570  if (x & C_FIXNUM_BIT)
571    C_kontinue(k, C_a_u_i_fixnum_abs(&a, 1, x));
572  else
573    C_u_bignum_abs(1, (C_word)NULL, k, x);
574}
575
576void C_ccall
577C_u_bignum_abs(C_word c, C_word self, C_word k, C_word x)
578{
579  if (C_bignum_negativep(x))
580    C_u_bignum_negate(1, (C_word)NULL, k, x);
581  else
582    C_kontinue(k, x);
583}
584
585void C_ccall
586C_basic_negate(C_word c, C_word self, C_word k, C_word x)
587{
588  C_word ab[nmax(C_SIZEOF_FIX_BIGNUM, C_SIZEOF_FLONUM)], *a = ab;
589
590  if (x & C_FIXNUM_BIT)
591    C_kontinue(k, C_a_u_i_fixnum_negate(&a, 1, x));
592  else if (C_immediatep(x))
593    C_kontinue(k, C_SCHEME_FALSE);
594  else if (C_block_header(x) == C_FLONUM_TAG)
595    C_kontinue(k, C_a_i_flonum_negate(&a, 1, x));
596  else if (C_IS_BIGNUM_TYPE(x))
597    C_u_bignum_negate(1, (C_word)NULL, k, x);
598  else
599    C_kontinue(k, C_SCHEME_FALSE);
600}
601
602void C_ccall
603C_u_integer_negate(C_word c, C_word self, C_word k, C_word x)
604{
605  C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
606
607  if (x & C_FIXNUM_BIT)
608    C_kontinue(k, C_a_u_i_fixnum_negate(&a, 1, x));
609  else
610    C_u_bignum_negate(1, (C_word)NULL, k, x);
611}
612
613C_regparm C_word C_fcall
614C_a_u_i_fixnum_negate(C_word **ptr, C_word n, C_word x)
615{
616  /* Exceptional situation: this will cause an overflow to itself */
617  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM)) /* C_fitsinfixnump(x) */
618    return C_bignum2(ptr, 0, 0, 1);
619  else
620    return C_fix(-C_unfix(x));
621}
622
623void C_ccall
624C_u_bignum_negate(C_word c, C_word self, C_word k, C_word x)
625{
626  if (C_bignum_negated_fitsinfixnump(x)) {
627    C_kontinue(k, C_fix(C_MOST_NEGATIVE_FIXNUM));
628  } else {
629    C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2,
630           negp = C_mk_bool(!C_bignum_negativep(x)),
631           size = C_fix(C_bignum_size(x));
632    k2 = C_closure(&ka, 3, (C_word)bignum_negate_2, k, x);
633    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
634  }
635}
636
637static void
638bignum_negate_2(C_word c, C_word self, C_word new_big)
639{
640  C_word k = C_block_item(self, 1),
641         old_big = C_block_item(self, 2);
642  bignum_digits_destructive_copy(new_big, old_big);
643  C_kontinue(k, C_bignum_normalize(new_big));
644}
645
646C_word
647C_u_i_bignum_cmp(C_word x, C_word y)
648{
649  if (C_bignum_negativep(x)) {
650    if (C_bignum_negativep(y)) { /* Largest negative number is smallest */
651      return C_fix(bignum_cmp_unsigned(y, x));
652    } else {
653      return C_fix(-1);
654    }
655  } else {
656    if (C_bignum_negativep(y)) {
657      return C_fix(1);
658    } else {
659      return C_fix(bignum_cmp_unsigned(x, y));
660    }
661  }
662}
663
664/* XXX TODO: Maybe pass true/false/bignum as initp, to allow copying data? */
665void C_ccall
666C_allocate_bignum(C_word c, C_word self, C_word k, C_word size, C_word negp, C_word initp)
667{
668  C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, init;
669  k2 = C_closure(&ka, 3, (C_word)allocate_bignum_2, k, negp);
670
671  init = C_and(initp, C_make_character('\0'));
672  C_allocate_vector(6, (C_word)NULL, k2,
673                    C_bytes(C_fixnum_plus(size, C_fix(1))), /* Add header */
674                    /* Byte vec, initialization, align at 8 bytes */
675                    C_SCHEME_TRUE, init, C_SCHEME_FALSE);
676}
677
678static void
679allocate_bignum_2(C_word c, C_word self, C_word bigvec)
680{
681  C_word ab[C_SIZEOF_STRUCTURE(2)], *a = ab, bignum,
682         k = C_block_item(self, 1),
683         negp = C_truep(C_block_item(self, 2)),
684         size = C_bytestowords(C_header_size(bigvec))-1;
685
686  C_word tagvec = CHICKEN_gc_root_ref(tags);
687
688  C_set_block_item(bigvec, 0, negp ? C_BIGNUM_HEADER_SIGN_BIT | size : size);
689
690  bignum = C_structure(&a, 2, C_block_item(tagvec, BIG_TAG), bigvec);
691  C_kontinue(k, bignum);
692}
693
694/* Normalization: scan trailing zeroes, then return a fixnum if the
695 * value fits, or trim the bignum's length. */
696C_word C_ccall
697C_bignum_normalize(C_word big)
698{
699  C_word *start = C_bignum_digits(big);
700  C_word *last_digit = start + C_bignum_size(big) - 1;
701  C_word *scan = last_digit, length;
702
703  while (scan >= start && *scan == 0)
704    scan--;
705  length = scan - start + 1;
706 
707  switch(length) {
708  case 0:
709    return C_fix(0);
710  case 1:
711    return C_fix(C_bignum_negativep(big) ? -*start : *start);
712  case 2:
713    if (C_bignum_negativep(big) && *scan == 1 && *start == 0)
714      return C_fix(C_MOST_NEGATIVE_FIXNUM);
715    /* FALLTHROUGH */
716  default:
717    if (scan < last_digit) {
718      /* Mutate vector size of internal bignum vector. */
719      C_block_header(C_internal_bignum(big)) = (C_STRING_TYPE | C_wordstobytes(length+1));
720      /* Set internal header. */
721      C_bignum_header(big) = (C_bignum_header(big) & C_BIGNUM_HEADER_SIGN_BIT) | length;
722    }
723    return big;
724  }
725}
726
727static C_word
728bignum_digits_destructive_scale_up_with_carry(C_word *start, C_word *end, C_word factor, C_word carry)
729{
730  C_word digit, product_hi, product_lo, *scan = start;
731
732  assert(C_fitsinbignumhalfdigitp(carry));
733
734  while (scan < end) {
735    digit = (*scan);
736    product_lo = factor * C_BIGNUM_DIGIT_LO_HALF(digit) + carry;
737    product_hi = factor * C_BIGNUM_DIGIT_HI_HALF(digit) +
738            C_BIGNUM_DIGIT_HI_HALF(product_lo);
739    (*scan++) = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(product_hi),
740                                       C_BIGNUM_DIGIT_LO_HALF(product_lo));
741    carry = C_BIGNUM_DIGIT_HI_HALF(product_hi);
742  }
743  return carry;
744}
745
746/* Given (denominator > 1), it is fairly easy to show that
747   quotient_high fits a bignum halfdigit, after which it is easy to
748   see that all digits fit a bignum full digit.
749   
750   This works because denominator is known to fit a halfdigit, which
751   means that the remainder modulo denominator will also fit a halfdigit. */
752static C_word
753bignum_digits_destructive_scale_down(C_word *start, C_word *end, C_word denominator)
754{
755  C_word numerator, remainder = 0, digit, quotient_high, *scan = end;
756
757  assert((denominator > 1) && C_fitsinbignumhalfdigitp(denominator));
758  while (start <= scan) {
759    digit = *scan;
760    numerator = C_BIGNUM_DIGIT_COMBINE(remainder, C_BIGNUM_DIGIT_HI_HALF(digit));
761    quotient_high = (numerator / denominator);
762    numerator = C_BIGNUM_DIGIT_COMBINE(numerator % denominator,
763                                       C_BIGNUM_DIGIT_LO_HALF(digit));
764    (*scan--) = C_BIGNUM_DIGIT_COMBINE(quotient_high, numerator / denominator);
765    remainder = numerator % denominator;
766  }
767  return remainder;
768}
769
770void C_ccall
771C_2_basic_times(C_word c, C_word self, C_word k, C_word x, C_word y)
772{
773  C_word ab[nmax(/*C_SIZEOF_FIX_BIGNUM, */ C_SIZEOF_FLONUM * 2,
774                 C_SIZEOF_FIX_BIGNUM * 2 + C_SIZEOF_BIGNUM(4))],
775         *a = ab;
776
777  if (x & C_FIXNUM_BIT) {
778    if (y & C_FIXNUM_BIT) {
779      C_kontinue(k, C_a_u_i_2_fixnum_times(&a, 2, x, y));
780    } else if (C_immediatep(y)) {
781      C_kontinue(k, C_SCHEME_FALSE);
782    } else if (C_block_header(y) == C_FLONUM_TAG) {
783      C_kontinue(k, C_flonum(&a, (double)C_unfix(x) * C_flonum_magnitude(y)));
784    } else if (C_IS_BIGNUM_TYPE(y)) {
785      C_u_2_integer_times(2, (C_word)NULL, k, x, y);
786    } else {
787      C_kontinue(k, C_SCHEME_FALSE);
788    }
789  } else if (C_immediatep(x)) {
790    C_kontinue(k, C_SCHEME_FALSE);
791  } else if (C_block_header(x) == C_FLONUM_TAG) {
792    if (y & C_FIXNUM_BIT) {
793      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) * (double)C_unfix(y)));
794    } else if (C_immediatep(y)) {
795      C_kontinue(k, C_SCHEME_FALSE);
796    } else if (C_block_header(y) == C_FLONUM_TAG) {
797      C_kontinue(k, C_a_i_flonum_times(&a, 2, x, y));
798    } else if (C_IS_BIGNUM_TYPE(y)) {
799      C_kontinue(k, C_a_i_flonum_times(&a, 2, x, C_a_u_i_big_to_flo(&a, 1, y)));
800    } else {
801      C_kontinue(k, C_SCHEME_FALSE);
802    }
803  } else if (C_IS_BIGNUM_TYPE(x)) {
804    if (y & C_FIXNUM_BIT) {
805      C_u_2_integer_times(2, (C_word)NULL, k, x, y);
806    } else if (C_immediatep(y)) {
807      C_kontinue(k, C_SCHEME_FALSE);
808    } else if (C_block_header(y) == C_FLONUM_TAG) {
809      C_kontinue(k, C_a_i_flonum_times(&a, 2, C_a_u_i_big_to_flo(&a, 1, x), y));
810    } else if (C_IS_BIGNUM_TYPE(y)) {
811      C_u_2_bignum_times(2, (C_word)NULL, k, x, y);
812    } else {
813      C_kontinue(k, C_SCHEME_FALSE);
814    }
815  } else {
816    C_kontinue(k, C_SCHEME_FALSE);
817  }
818}
819
820/* XXX TODO: Rework to allocate C_SIZEOF_BIGNUM(4)?  Not worth the hassle? */
821/* Needs at most 2 * SIZEOF_FIX_BIGNUM + SIZEOF_BIGNUM(4) */
822C_regparm C_word C_fcall
823C_a_u_i_2_fixnum_times(C_word **ptr, C_word n, C_word x, C_word y)
824{
825  C_word absx, absy, negp, *d, r;
826
827  /* We don't strictly need the abses in all branches... */
828  absx = C_unfix(x);
829  absx = absx < 0 ? -absx : absx;
830  absy = C_unfix(y);
831  absy = absy < 0 ? -absy : absy;
832  negp = ((x & C_INT_SIGN_BIT) ? !(y & C_INT_SIGN_BIT) : (y & C_INT_SIGN_BIT));
833
834  if (C_fitsinbignumhalfdigitp(absx)) {
835     if (x == C_fix(0) || x == C_fix(1) || C_fitsinbignumhalfdigitp(absy)) {
836       return C_fix(negp ? -(absx * absy) : (absx * absy));
837     } else {
838       if (y == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
839         y = C_bignum2(ptr, negp != 0, 0, 1); /* Two is always enough */
840       } else {
841         y = C_bignum2(ptr, negp != 0, absy, 0); /* May need one for carry */
842       }
843       d = C_bignum_digits(y);
844       r = bignum_digits_destructive_scale_up_with_carry(d, d+2, absx, 0);
845       assert(r == 0); /* Should never result in a carry; y is big enough */
846       return C_bignum_normalize(y);
847     }
848  } else if (C_fitsinbignumhalfdigitp(absy)) {
849     if (absy == 0 || y == C_fix(1) /*|| C_fitsinbignumhalfdigitp(absx) */) {
850       return C_fix(negp ? -(absx * absy) : (absx * absy));
851     } else {
852       if (x == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
853         x = C_bignum2(ptr, negp != 0, 0, 1); /* Two is always enough */
854       } else {
855         x = C_bignum2(ptr, negp != 0, absx, 0); /* May need one for carry */
856       }
857       d = C_bignum_digits(x);
858       r = bignum_digits_destructive_scale_up_with_carry(d, d+2, absy, 0);
859       assert(r == 0); /* Should never result in a carry; x is big enough */
860       return C_bignum_normalize(x);
861     }
862  } else {
863    x = C_a_u_i_fix_to_big(ptr, x);
864    y = C_a_u_i_fix_to_big(ptr, y);
865    r = C_bignum4(ptr, negp != 0, 0, 0, 0, 0);
866    bignum_digits_multiply(x, y, r);
867    return C_bignum_normalize(r);
868  }
869}
870
871void C_ccall
872C_u_2_integer_times(C_word c, C_word self, C_word k, C_word x, C_word y)
873{
874  C_word ab[nmax(/* C_SIZEOF_FIX_BIGNUM, */ C_SIZEOF_CLOSURE(4),
875                 C_SIZEOF_FIX_BIGNUM * 2 + C_SIZEOF_BIGNUM(4))],
876         *a = ab, absy, negp, size, k2;
877
878  if (x & C_FIXNUM_BIT) {
879    if (y & C_FIXNUM_BIT) {
880      C_kontinue(k, C_a_u_i_2_fixnum_times(&a, 2, x, y));
881    } else {
882      absy = x; /* swap to ensure x is a bignum and y a fixnum */
883      x = y;
884      y = absy;
885    }
886  }
887  /* Here, we know for sure that X is a bignum */
888  if (y & C_FIXNUM_BIT) {
889    absy = C_unfix(y);
890    absy = absy < 0 ? -absy : absy;
891    negp = C_mk_bool((y & C_INT_SIGN_BIT) ?
892                     !C_bignum_negativep(x) :
893                     C_bignum_negativep(x));
894 
895    if (C_fitsinbignumhalfdigitp(absy)) {
896      k2 = C_closure(&a, 4, (C_word)integer_times_2, k, x, C_fix(absy));
897      size = C_fix(C_bignum_size(x) + 1); /* Needs _at most_ one more digit */
898      C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
899    } else {
900      y = C_a_u_i_fix_to_big(&a, y);
901      bignum_times_bignum_unsigned(k, x, y, negp);
902    }
903  } else {
904    C_u_2_bignum_times(2, (C_word)NULL, k, x, y);
905  }
906}
907
908static void
909integer_times_2(C_word c, C_word self, C_word new_big)
910{
911  C_word k = C_block_item(self, 1),
912         old_bigx = C_block_item(self, 2),
913         fixy = C_block_item(self, 3),
914         *digits = C_bignum_digits(new_big),
915         *end_digit = digits + C_bignum_size(old_bigx);
916
917  bignum_digits_destructive_copy(new_big, old_bigx);
918
919  /* Scale up, and sanitise the result. */
920  *end_digit = bignum_digits_destructive_scale_up_with_carry(digits, end_digit,
921                                                             C_unfix(fixy), 0);
922  C_kontinue(k, C_bignum_normalize(new_big));
923}
924
925void C_ccall
926C_u_2_bignum_times(C_word c, C_word self, C_word k, C_word x, C_word y)
927{
928  C_word negp = C_bignum_negativep(x) ?
929                !C_bignum_negativep(y) :
930                C_bignum_negativep(y);
931  bignum_times_bignum_unsigned(k, x, y, C_mk_bool(negp));
932}
933
934/* Multiplication
935   Maximum value for product_lo or product_hi:
936        ((R * R) + (R * (R - 2)) + (R - 1))
937   Maximum value for carry: ((R * (R - 1)) + (R - 1))
938        where R == 2^HALF_DIGIT_LENGTH */
939static void
940bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp)
941{
942  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
943 
944  if (C_bignum_size(y) > C_bignum_size(x)) { /* Ensure size(x) <= size(y) */
945    C_word z = x;
946    x = y;
947    y = z;
948  }
949
950  k2 = C_closure(&ka, 4, (C_word)bignum_times_bignum_unsigned_2, k, x, y);
951
952  size = C_fix(C_bignum_size(x) + C_bignum_size(y));
953  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
954}
955
956static C_regparm void
957bignum_digits_multiply(C_word x, C_word y, C_word result)
958{
959  C_word carry, y_digit_lo, y_digit_hi, x_digit_lo,
960         x_digit_hi, product_lo, *scan_r, *scan_y,
961         x_digit, y_digit, product_hi,
962         *scan_x = C_bignum_digits(x),
963         *end_x = scan_x + C_bignum_size(x),
964         *start_y = C_bignum_digits(y),
965         *end_y = start_y + C_bignum_size(y),
966         *start_r = C_bignum_digits(result);
967
968  while (scan_x < end_x) {
969    x_digit = (*scan_x++);
970    x_digit_lo = C_BIGNUM_DIGIT_LO_HALF(x_digit);
971    x_digit_hi = C_BIGNUM_DIGIT_HI_HALF(x_digit);
972    carry = 0;
973    scan_y = start_y;
974    scan_r = (start_r++);
975
976    while (scan_y < end_y) {
977      y_digit = (*scan_y++);
978      y_digit_lo = C_BIGNUM_DIGIT_LO_HALF(y_digit);
979      y_digit_hi = C_BIGNUM_DIGIT_HI_HALF(y_digit);
980
981      product_lo = (*scan_r) +
982                   x_digit_lo * y_digit_lo +
983                   C_BIGNUM_DIGIT_LO_HALF(carry);
984
985      product_hi = x_digit_hi * y_digit_lo +
986                   x_digit_lo * y_digit_hi +
987                   C_BIGNUM_DIGIT_HI_HALF(product_lo) +
988                   C_BIGNUM_DIGIT_HI_HALF(carry);
989
990      (*scan_r++) = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(product_hi),
991                                           C_BIGNUM_DIGIT_LO_HALF(product_lo));
992
993      carry = x_digit_hi * y_digit_hi + C_BIGNUM_DIGIT_HI_HALF(product_hi);
994    }
995    (*scan_r) += carry;
996  }
997}
998
999static void
1000bignum_times_bignum_unsigned_2(C_word c, C_word self, C_word result)
1001{
1002  C_word k = C_block_item(self, 1),
1003         x = C_block_item(self, 2),
1004         y = C_block_item(self, 3);
1005
1006  bignum_digits_multiply(x, y, result);
1007  C_kontinue(k, C_bignum_normalize(result));
1008}
1009
1010void C_ccall
1011C_digits_to_integer(C_word c, C_word self, C_word k, C_word str,
1012                    C_word start, C_word end, C_word radix, C_word negp)
1013{
1014  C_word kab[C_SIZEOF_CLOSURE(6)], *ka = kab, k2, size;
1015  size_t nbits;
1016
1017  assert((C_unfix(radix) > 1) && C_fitsinbignumhalfdigitp(C_unfix(radix)));
1018 
1019  if (start == end) {
1020    C_kontinue(k, C_SCHEME_FALSE);
1021  } else {
1022    k2 = C_closure(&ka, 6, (C_word)digits_to_integer_2, k, str, start, end, radix);
1023 
1024    nbits = (C_unfix(end) - C_unfix(start)) * ilen(C_unfix(radix));
1025    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(nbits));
1026    /* XXX: Why initialize? */
1027    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
1028  }
1029}
1030
1031/* XXX TODO: This should be faster, dammit! */
1032static void
1033digits_to_integer_2(C_word c, C_word self, C_word result)
1034{
1035  C_word k = C_block_item(self, 1),
1036         str = C_block_item(self, 2),
1037         start = C_unfix(C_block_item(self, 3)),
1038         end = C_unfix(C_block_item(self, 4)),
1039         radix = C_unfix(C_block_item(self, 5)),
1040         first_digit = C_unfix(C_block_item(self, 6)),
1041         *digits = C_bignum_digits(result),
1042         *last_digit = digits, /* Initially, bignum is all zeroes */
1043         big_digit = 0, factor = radix,
1044         next_big_digit, next_factor,
1045         carry, str_digit;
1046  char *str_scan = C_c_string(str) + start,
1047       *str_end = C_c_string(str) + end;
1048
1049  /* Hash characters in numbers are mapped to 0 */
1050# define HEXDIGIT_CHAR_TO_INT(x)                                        \
1051    (((x) == '#') ? 0 :                                                 \
1052     (((x) >= (int)'a') ? ((x) - (int)'a' + 10) : ((x) - (int)'0')))
1053
1054  /* This tries to save up as much as possible in the local C_word
1055   * big_digit, and only when it exceeds what we would be able to
1056   * multiply easily, we scale up the bignum and add what we saved up.
1057   */
1058  while (str_scan < str_end) {
1059    str_digit = HEXDIGIT_CHAR_TO_INT(C_tolower((int)*str_scan));
1060    str_scan++;
1061
1062    next_big_digit = big_digit * radix;
1063    next_big_digit += str_digit;
1064    next_factor = factor * radix;
1065
1066    if (str_digit >= radix || str_digit < 0) {
1067      C_kontinue(k, C_SCHEME_FALSE);
1068    } else if (C_fitsinbignumhalfdigitp(next_big_digit) &&
1069               C_fitsinbignumhalfdigitp(next_factor)) {
1070      factor = next_factor;
1071      big_digit = next_big_digit;
1072    } else {
1073      carry = bignum_digits_destructive_scale_up_with_carry(
1074              digits, last_digit, factor, big_digit);
1075
1076      if (carry) (*last_digit++) = carry; /* Move end */
1077
1078      big_digit = str_digit;
1079      factor = radix;
1080    }
1081  }
1082# undef HEXDIGIT_CHAR_TO_INT
1083
1084  /* Final step.  We always must do this, because the loop never
1085   * processes the "current" character into the bignum (lookahead 1).
1086   */
1087  carry = bignum_digits_destructive_scale_up_with_carry(
1088          digits, last_digit, factor, big_digit);
1089  if (carry) (*last_digit++) = carry; /* Move end */
1090
1091  C_kontinue(k, C_bignum_normalize(result));
1092}
1093
1094void C_ccall
1095C_u_bignum_to_digits(C_word c, C_word self, C_word k, C_word value, C_word radix)
1096{
1097  /* This copies the bignum over into a working copy that can be mutated. */
1098  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size,
1099         negp = C_mk_bool(C_bignum_negativep(value));
1100
1101  k2 = C_closure(&ka, 4, (C_word)bignum_to_digits_2, k, value, radix);
1102
1103  size = C_fix(C_bignum_size(value));
1104  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1105}
1106
1107static void
1108bignum_to_digits_2(C_word c, C_word self, C_word working_copy)
1109{
1110  C_word k = C_block_item(self, 1),
1111         old_big = C_block_item(self, 2),
1112         radix = C_unfix(C_block_item(self, 3)),
1113         len;
1114
1115  bignum_digits_destructive_copy(working_copy, old_big);
1116
1117  assert(radix > 1 && C_fitsinbignumhalfdigitp(radix));
1118
1119  /* A sloppy over-approximation of the number of radix digits we
1120   * need.  Due to the algorithm chopping off slices of "base" size at
1121   * a time, we add one to the bignum size, so we are sure we can
1122   * handle any leading zeroes.  The extra 1 we add at the end is for
1123   * an optional sign.
1124   */
1125  len = (C_bignum_size(old_big)+1)*C_BIGNUM_DIGIT_LENGTH / (ilen(radix)-1) + 1;
1126
1127  /* Nice: We can recycle the current closure */
1128  C_set_block_item(self, 0, (C_word)bignum_to_digits_3);
1129  /* item 1 is still the same continuation */
1130  C_set_block_item(self, 2, working_copy);
1131
1132  C_allocate_vector(6, (C_word)NULL, self, C_fix(len),
1133                    /* Byte vec, no initialization, align at 8 bytes */
1134                    C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
1135}
1136
1137/* XXX TODO: This should be MUCH faster, dammit! */
1138static void
1139bignum_to_digits_3(C_word c, C_word self, C_word string)
1140{
1141  static char *characters = "0123456789abcdef";
1142  C_word k = C_block_item(self, 1),
1143         working_copy = C_block_item(self, 2),
1144         radix = C_unfix(C_block_item(self, 3)),
1145         *start = C_bignum_digits(working_copy),
1146         *scan = start + C_bignum_size(working_copy) - 1,
1147         len = C_header_size(string),
1148         digit, steps, i, base;
1149  char *buf = C_c_string(string), *index = buf + len - 1;
1150
1151  /* Calculate the largest power of radix that fits a halfdigit:
1152   * steps = log10(2^halfdigit_bits), base = 10^steps
1153   */
1154  for(steps = 0, base = radix; C_fitsinbignumhalfdigitp(base); base *= radix)
1155    steps++;
1156
1157  base /= radix; /* Back down: we overshot in the loop */
1158
1159  while (start <= scan) {
1160    digit = bignum_digits_destructive_scale_down(start, scan, base);
1161
1162    if (*scan == 0) scan--; /* Adjust if we exhausted the highest digit */
1163
1164    for(i = 0; i < steps; ++i) {
1165      *index-- = characters[digit % radix];
1166      digit /= radix;
1167    }
1168  }
1169  assert(index >= buf-1);
1170
1171  /* Move index onto first nonzero digit (advance, and skip leading zeroes) */
1172  while(*++index == '0');
1173
1174  if (C_bignum_negativep(working_copy))
1175    *--index = '-';
1176
1177  len -= index - buf; /* Shorten with distance between start and index. */
1178  C_memmove(buf, index, len); /* Move start of number to begin of string. */
1179  C_block_header(string) = C_STRING_TYPE | len; /* Mutate string length. */
1180  C_kontinue(k, string);
1181}
1182
1183C_word
1184C_a_u_i_big_to_flo(C_word **p, C_word n, C_word bignum)
1185{
1186  double accumulator = 0;
1187  C_word *start = C_bignum_digits(bignum);
1188  C_word *scan = start + C_bignum_size(bignum);
1189  while (start < scan) {
1190    accumulator *= (C_word)1 << C_BIGNUM_DIGIT_LENGTH;
1191    accumulator += (*--scan);
1192  }
1193  return C_flonum(p, (C_bignum_negativep(bignum) ? -accumulator : accumulator));
1194}
1195
1196void C_ccall
1197C_u_flo_to_int(C_word c, C_word self, C_word k, C_word x)
1198{
1199  int exponent;
1200  double significand = frexp(C_flonum_magnitude(x), &exponent);
1201
1202  if (exponent <= 0) {
1203    C_kontinue(k, C_fix(0));
1204  } else if (exponent == 1) { /* TODO: check significand * 2^exp fits fixnum? */
1205    C_kontinue(k, significand < 0.0 ? C_fix(-1) : C_fix(1));
1206  } else {
1207    C_word kab[C_SIZEOF_CLOSURE(4) + C_SIZEOF_FLONUM], *ka = kab, k2, size,
1208           negp = C_mk_bool(C_flonum_magnitude(x) < 0.0),
1209           sign = C_flonum(&ka, fabs(significand));
1210
1211    k2 = C_closure(&ka, 4, (C_word)flo_to_int_2, k, C_fix(exponent), sign);
1212
1213    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(exponent));
1214    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1215  }
1216}
1217
1218static void
1219flo_to_int_2(C_word c, C_word self, C_word result)
1220{
1221  C_word digit, k = C_block_item(self, 1),
1222         exponent = C_unfix(C_block_item(self, 2)),
1223         odd_bits = exponent % C_BIGNUM_DIGIT_LENGTH,
1224         *start = C_bignum_digits(result),
1225         *scan = start + C_bignum_size(result);
1226  /* It's always true that 0.5 <= s < 1 */
1227  double significand = C_flonum_magnitude(C_block_item(self, 3));
1228
1229  if (odd_bits > 0) { /* Handle most significant digit first */
1230    significand *= (C_word)1 << odd_bits;
1231    digit = (C_word)significand;
1232    (*--scan) = digit;
1233    significand -= (double)digit;
1234  }
1235
1236  while (start < scan && significand > 0) {
1237    significand *= (C_word)1 << C_BIGNUM_DIGIT_LENGTH;
1238    digit = (C_word)significand;
1239    (*--scan) = digit;
1240    significand -= (double)digit;
1241  }
1242
1243  /* Finish up by clearing any remaining, lower, digits */
1244  while (start < scan)
1245    (*--scan) = 0;
1246
1247  C_kontinue(k, C_bignum_normalize(result));
1248}
1249
1250/*
1251 * From Hacker's Delight by Henry S. Warren
1252 * based on a modified nlz() from section 5-3 (fig. 5-7)
1253 */
1254static C_word
1255ilen(C_uword x)
1256{
1257  C_uword y;
1258  C_word n = 0;
1259
1260#ifdef C_SIXTY_FOUR
1261  y = x >> 32; if (y != 0) { n += 32; x = y; }
1262#endif
1263  y = x >> 16; if (y != 0) { n += 16; x = y; }
1264  y = x >>  8; if (y != 0) { n +=  8; x = y; }
1265  y = x >>  4; if (y != 0) { n +=  4; x = y; }
1266  y = x >>  2; if (y != 0) { n +=  2; x = y; }
1267  y = x >>  1; if (y != 0) return n + 2;
1268  return n + x;
1269}
1270
1271C_word
1272C_u_i_int_length(C_word x)
1273{
1274  if (x & C_FIXNUM_BIT) {
1275    x = C_unfix(x);
1276    return C_fix(ilen((x < 0) ? ~x : x));
1277  } else { /* bignum */
1278    C_word len_1 = C_bignum_size(x) - 1,
1279           result = len_1 * C_BIGNUM_DIGIT_LENGTH,
1280           *startx = C_bignum_digits(x),
1281           *last_digit = C_bignum_digits(x) + len_1,
1282           last_digit_length = ilen(*last_digit);
1283
1284    /* If *only* the highest bit is set, negating results in one less bit */
1285    if (C_bignum_negativep(x) && *last_digit == (1 << (last_digit_length-1))) {
1286      while(startx < last_digit && *startx == 0) ++startx;
1287      if (startx == last_digit) result--;
1288    }
1289    return C_fix(result + last_digit_length);
1290  }
1291}
1292
1293void C_ccall
1294C_u_int_shift_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
1295{
1296  C_word kab[C_SIZEOF_FIX_BIGNUM * 2 + C_SIZEOF_CLOSURE(3) + C_SIZEOF_CLOSURE(2)],
1297         *ka = kab, k2, k3, size, minus1;
1298
1299  if (y == C_fix(0) || x == C_fix(0)) { /* Done (no shift) */
1300    C_kontinue(k, x);
1301  } else if (x & C_FIXNUM_BIT) {
1302    /* TODO: This should probably see if shifting fits a fixnum, first */
1303    x = C_a_u_i_fix_to_big(&ka, x);
1304  }
1305
1306  /* Invert all the bits before shifting right a negative value */
1307  if (C_bignum_negativep(x) && C_unfix(y) < 0) {
1308    /* When done shifting, invert again */
1309    k3 = C_closure(&ka, 2, (C_word)bignum_negate_after_shift, k);
1310    /* Before shifting, allocate the bignum */
1311    k2 = C_closure(&ka, 3, (C_word)bignum_allocate_for_shift, k3, y);
1312    /* Actually invert by subtracting: -1 - x */
1313    minus1 = C_a_u_i_fix_to_big(&ka, C_fix(-1));
1314    C_u_2_bignum_minus(2, (C_word)NULL, k2, minus1, x);
1315  } else {
1316    k2 = C_closure(&ka, 3, (C_word)bignum_allocate_for_shift, k, y);
1317    C_kontinue(k2, x);
1318  }
1319}
1320
1321static void
1322bignum_allocate_for_shift(C_word c, C_word self, C_word x)
1323{
1324  C_word k = C_block_item(self, 1),
1325         y = C_block_item(self, 2),
1326         uy = C_unfix(y),
1327         negp, digit_offset, bit_offset,
1328         ab[C_SIZEOF_FIX_BIGNUM + C_SIZEOF_CLOSURE(6)], *a = ab, k2, size;
1329
1330  if (x & C_FIXNUM_BIT) /* Normalisation may happen after negation */
1331    x = C_a_u_i_fix_to_big(&a, x);
1332
1333  negp = C_mk_bool(C_bignum_negativep(x));
1334 
1335  /* uy is guaranteed not to be 0 here */
1336  if (uy > 0) {
1337    digit_offset = uy / C_BIGNUM_DIGIT_LENGTH;
1338    bit_offset =   uy % C_BIGNUM_DIGIT_LENGTH;
1339
1340    k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k,
1341                   x, C_SCHEME_TRUE, C_fix(digit_offset), C_fix(bit_offset));
1342    size = C_fix(C_bignum_size(x) + digit_offset + 1);
1343    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
1344  } else if (-uy >= C_bignum_size(x) * (C_word)C_BIGNUM_DIGIT_LENGTH) {
1345    /* All bits are shifted out, just return 0 */
1346    C_kontinue(k, C_fix(0));
1347  } else {
1348    digit_offset = -uy / C_BIGNUM_DIGIT_LENGTH;
1349    bit_offset =   -uy % C_BIGNUM_DIGIT_LENGTH;
1350   
1351    k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k,
1352                   x, C_SCHEME_FALSE, C_fix(digit_offset), C_fix(bit_offset));
1353    size = C_fix(C_bignum_size(x) - digit_offset);
1354    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
1355  }
1356}
1357
1358static void
1359bignum_negate_after_shift(C_word c, C_word self, C_word result)
1360{
1361  C_word k = C_block_item(self, 1);
1362  if (result & C_FIXNUM_BIT) { /* Normalisation may happen after shift */
1363    C_kontinue(k, C_fix(-1 - C_unfix(result)));
1364  } else {
1365    C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab, minus1;
1366    minus1 = C_a_u_i_fix_to_big(&a, C_fix(-1));
1367    C_u_2_bignum_minus(2, (C_word)NULL, k, minus1, result);
1368  }
1369}
1370
1371static void
1372bignum_actual_shift(C_word c, C_word self, C_word result)
1373{
1374  C_word k = C_block_item(self, 1),
1375         x = C_block_item(self, 2),
1376         shift_left = C_truep(C_block_item(self, 3)),
1377         digit_offset = C_unfix(C_block_item(self, 4)),
1378         bit_offset = C_unfix(C_block_item(self, 5)),
1379         *scanx, *scanr, *end;
1380
1381  if (shift_left) {
1382    scanr = C_bignum_digits(result) + digit_offset;
1383    scanx = C_bignum_digits(x);
1384    end = scanx + C_bignum_size(x);
1385   
1386    while (scanx < end) {
1387      *scanr = *scanr | (*scanx & C_BIGNUM_DIGIT_MASK) << bit_offset;
1388      *scanr = *scanr & C_BIGNUM_DIGIT_MASK;
1389      scanr++;
1390      *scanr = *scanx++ >> (C_BIGNUM_DIGIT_LENGTH - bit_offset);
1391      *scanr = *scanr & C_BIGNUM_DIGIT_MASK;
1392    }
1393  } else {
1394    scanr = C_bignum_digits(result);
1395    scanx = C_bignum_digits(x) + digit_offset;
1396    end = scanr + C_bignum_size(result) - 1;
1397   
1398    while (scanr < end) {
1399      *scanr =  (*scanx++ & C_BIGNUM_DIGIT_MASK) >> bit_offset;
1400      *scanr = (*scanr | 
1401        *scanx << (C_BIGNUM_DIGIT_LENGTH - bit_offset)) & C_BIGNUM_DIGIT_MASK;
1402      scanr++;
1403    }
1404    *scanr =  (*scanx++ & C_BIGNUM_DIGIT_MASK) >> bit_offset;
1405  }
1406  C_kontinue(k, C_bignum_normalize(result));
1407}
1408
1409/*
1410 * This random number generator is very simple. Probably too simple...
1411 */
1412C_word C_ccall
1413C_u_i_bignum_randomize(C_word bignum)
1414{
1415  C_word seed = 0,
1416         *scan = C_bignum_digits(bignum),
1417         *end = scan + C_bignum_size(bignum);
1418
1419  /* What a cheap way to initialize the random generator. I feel dirty! */
1420  while (scan < end)
1421    seed ^= *scan++;
1422
1423  srand(seed);
1424  return C_SCHEME_UNDEFINED;
1425}
1426
1427void C_ccall
1428C_u_bignum_random(C_word c, C_word self, C_word k, C_word max)
1429{
1430  C_word k2, kab[C_SIZEOF_CLOSURE(4)], *ka = kab, size,
1431         max_len, max_bits, max_top_digit, d, negp;
1432
1433  max_len = C_bignum_size(max);
1434  max_top_digit = d = C_bignum_digits(max)[max_len - 1];
1435 
1436  max_bits = (max_len - 1) * C_BIGNUM_DIGIT_LENGTH;
1437  while(d) {
1438    max_bits++;
1439    d >>= 1;
1440  }
1441  /* Subtract/add one because we don't want zero to be over-represented */
1442  size = ((double)rand())/(RAND_MAX + 1.0) * (double)(max_bits - 1);
1443  size = C_fix(C_BIGNUM_BITS_TO_DIGITS(size + 1));
1444
1445  negp = C_mk_bool(C_bignum_negativep(max));
1446  k2 = C_closure(&ka, 4, (C_word)bignum_random_2, k, C_fix(max_top_digit), C_fix(max_len));
1447  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1448}
1449
1450static void
1451bignum_random_2(C_word c, C_word self, C_word result)
1452{
1453  C_word k = C_block_item(self, 1),
1454         max_top_digit = C_unfix(C_block_item(self, 2)),
1455         max_len = C_unfix(C_block_item(self, 3)),
1456         *scan = C_bignum_digits(result),
1457         *end = scan + C_bignum_size(result); /* Go to just before the end. */
1458
1459  while(scan < end)
1460    *scan++ = ((double)rand())/(RAND_MAX + 1.0) * (double)((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
1461  /*
1462   * Last word is special when length is max_len: It must be less than
1463   * max's most significant digit, instead of BIGNUM_RADIX.
1464   */
1465  if (max_len == C_bignum_size(result))
1466    *scan = ((double)rand())/(RAND_MAX + 1.0) * (double)max_top_digit;
1467  else
1468    *scan = ((double)rand())/(RAND_MAX + 1.0) * (double)((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
1469
1470  C_kontinue(k, C_bignum_normalize(result));
1471}
1472
1473
1474/* op identifies the operator: 0 means AND, 1 means IOR, 2 means XOR */
1475void C_ccall
1476C_u_2_integer_bitwise_op(C_word c, C_word self, C_word k, C_word op, C_word x, C_word y)
1477{
1478  C_word kab[C_SIZEOF_FIX_BIGNUM*2 + C_SIZEOF_CLOSURE(5)], *ka = kab, k2,
1479         size, negp;
1480 
1481  if (x & C_FIXNUM_BIT)
1482    x = C_a_u_i_fix_to_big(&ka, x);
1483
1484  if (y & C_FIXNUM_BIT)
1485    y = C_a_u_i_fix_to_big(&ka, y);
1486
1487  if (C_bignum_negativep(x)) {
1488    if (C_bignum_negativep(y)) {
1489      negp = C_mk_bool(op == C_fix(0) || op == C_fix(1)); /* and / ior */
1490      size = C_fix(nmax(C_bignum_size(x) + 1, C_bignum_size(y) + 1));
1491      k2 = C_closure(&ka, 5, (C_word)bignum_negneg_bitwise_op, k, op, x, y);
1492    } else {
1493      negp = C_mk_bool(op == C_fix(1) || op == C_fix(2)); /* ior / xor */
1494      size = C_fix(nmax(C_bignum_size(y), C_bignum_size(x) + 1)); /*!*/
1495      k2 = C_closure(&ka, 5, (C_word)bignum_posneg_bitwise_op, k, op, y, x);
1496    }
1497  } else {
1498    if (C_bignum_negativep(y)) {
1499      negp = C_mk_bool(op == C_fix(1) || op == C_fix(2)); /* ior / xor */
1500      size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y) + 1));
1501      k2 = C_closure(&ka, 5, (C_word)bignum_posneg_bitwise_op, k, op, x, y);
1502    } else {
1503      negp = C_SCHEME_FALSE;
1504      size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y)));
1505      k2 = C_closure(&ka, 5, (C_word)bignum_pospos_bitwise_op, k, op, x, y);
1506    }
1507  }
1508
1509  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1510}
1511
1512static void
1513bignum_pospos_bitwise_op(C_word c, C_word self, C_word result)
1514{
1515  C_word k = C_block_item(self, 1),
1516         op = C_block_item(self, 2),
1517         arg1 = C_block_item(self, 3),
1518         arg2 = C_block_item(self, 4),
1519         *scanr = C_bignum_digits(result),
1520         *endr = scanr + C_bignum_size(result),
1521         *scan1 = C_bignum_digits(arg1),
1522         *end1 = scan1 + C_bignum_size(arg1),
1523         *scan2 = C_bignum_digits(arg2),
1524         *end2 = scan2 + C_bignum_size(arg2),
1525         digit1, digit2;
1526
1527  while (scanr < endr) {
1528    digit1 = (scan1 < end1) ? *scan1++ : 0;
1529    digit2 = (scan2 < end2) ? *scan2++ : 0;
1530    *scanr++ = (op == C_fix(0)) ? digit1 & digit2 :
1531               (op == C_fix(1)) ? digit1 | digit2 :
1532                                  digit1 ^ digit2;
1533  }
1534  C_kontinue(k, C_bignum_normalize(result));
1535}
1536
1537static void
1538bignum_posneg_bitwise_op(C_word c, C_word self, C_word result)
1539{
1540  C_word k = C_block_item(self, 1),
1541         op = C_block_item(self, 2),
1542         arg1 = C_block_item(self, 3),
1543         arg2 = C_block_item(self, 4),
1544         *scanr = C_bignum_digits(result),
1545         *endr = scanr + C_bignum_size(result),
1546         *scan1 = C_bignum_digits(arg1),
1547         *end1 = scan1 + C_bignum_size(arg1),
1548         *scan2 = C_bignum_digits(arg2),
1549         *end2 = scan2 + C_bignum_size(arg2),
1550         digit1, digit2, carry2 = 1;
1551
1552  while (scanr < endr) {
1553    digit1 = (scan1 < end1) ? *scan1++ : 0;
1554    digit2 = (~((scan2 < end2) ? *scan2++ : 0) & C_BIGNUM_DIGIT_MASK) + carry2;
1555
1556    if (C_fitsinbignumdigitp(digit2)) {
1557      carry2 = 0;
1558    } else {
1559      digit2 &= C_BIGNUM_DIGIT_MASK;
1560      carry2 = 1;
1561    }
1562   
1563    *scanr++ = (op == C_fix(0)) ? digit1 & digit2 :
1564               (op == C_fix(1)) ? digit1 | digit2 :
1565                                  digit1 ^ digit2;
1566  }
1567
1568  bignum_maybe_negate_magnitude(k, result);
1569}
1570
1571static void
1572bignum_negneg_bitwise_op(C_word c, C_word self, C_word result)
1573{
1574  C_word k = C_block_item(self, 1),
1575         op = C_block_item(self, 2),
1576         arg1 = C_block_item(self, 3),
1577         arg2 = C_block_item(self, 4),
1578         *scanr = C_bignum_digits(result),
1579         *endr = scanr + C_bignum_size(result),
1580         *scan1 = C_bignum_digits(arg1),
1581         *end1 = scan1 + C_bignum_size(arg1),
1582         *scan2 = C_bignum_digits(arg2),
1583         *end2 = scan2 + C_bignum_size(arg2),
1584         digit1, digit2, carry1 = 1, carry2 = 1;
1585
1586  while (scanr < endr) {
1587    digit1 = (~((scan1 < end1) ? *scan1++ : 0) & C_BIGNUM_DIGIT_MASK) + carry1;
1588    digit2 = (~((scan2 < end2) ? *scan2++ : 0) & C_BIGNUM_DIGIT_MASK) + carry2;
1589
1590    if (C_fitsinbignumdigitp(digit1)) {
1591      carry1 = 0;
1592    } else {
1593      digit1 &= C_BIGNUM_DIGIT_MASK;
1594      carry1 = 1;
1595    }
1596    if (C_fitsinbignumdigitp(digit2)) {
1597      carry2 = 0;
1598    } else {
1599      digit2 &= C_BIGNUM_DIGIT_MASK;
1600      carry2 = 1;
1601    }
1602   
1603    *scanr++ = (op == C_fix(0)) ? digit1 & digit2 :
1604               (op == C_fix(1)) ? digit1 | digit2 :
1605                                  digit1 ^ digit2;
1606  }
1607
1608  bignum_maybe_negate_magnitude(k, result);
1609}
1610
1611static void
1612bignum_maybe_negate_magnitude(C_word k, C_word result)
1613{
1614  if (C_bignum_negativep(result)) {
1615    C_word *scan, *end, digit, carry;
1616
1617    scan = C_bignum_digits(result);
1618    end = scan + C_bignum_size(result);
1619    carry = 1;
1620
1621    while (scan < end) {
1622      digit = (~*scan & C_BIGNUM_DIGIT_MASK) + carry;
1623
1624      if (C_fitsinbignumdigitp(digit)) {
1625        carry = 0;
1626      } else {
1627        digit &= C_BIGNUM_DIGIT_MASK;
1628        carry = 1;
1629      }
1630   
1631      *scan++ = digit;
1632    }
1633  }
1634  C_kontinue(k, C_bignum_normalize(result));
1635}
1636
1637void C_ccall
1638C_u_bignum_divrem_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y)
1639{
1640  y = C_unfix(y);
1641
1642  if (y == 1) {
1643    C_values(4, C_SCHEME_UNDEFINED, k, x, C_fix(0));
1644  } else if (y == -1) {
1645    C_word kab[C_SIZEOF_CLOSURE(2)], *ka = kab, k2;
1646    k2 = C_closure(&ka, 2, (C_word)bignum_divrem_fixnum_2, k);
1647    C_u_bignum_negate(1, (C_word)NULL, k2, x);
1648  } else if (y == C_MOST_NEGATIVE_FIXNUM) {
1649    /* This is the only case we need to go allocate a bignum for */
1650    C_word kab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *ka = kab, k2, size,
1651           q_negp = C_mk_bool(!(C_bignum_negativep(x))), /* We know y < 0 */
1652           r_negp = C_mk_bool(C_bignum_negativep(x));
1653
1654    y = C_a_u_i_fix_to_big(&ka, C_fix(C_MOST_NEGATIVE_FIXNUM));
1655
1656    k2 = C_closure(&ka, 9, (C_word)bignum_divide_2_unsigned, k, x, y,
1657                   /* Return quotient *and* remainder */
1658                   C_SCHEME_TRUE, C_SCHEME_TRUE, r_negp,
1659                   /* Will be filled in later */
1660                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
1661    size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
1662    C_allocate_bignum(3, (C_word)NULL, k2, size, q_negp, C_SCHEME_FALSE);
1663  } else {
1664    C_word q_negp = C_mk_bool((y < 0) ?
1665                              !(C_bignum_negativep(x)) :
1666                              C_bignum_negativep(x)),
1667           r_negp = C_mk_bool(C_bignum_negativep(x)),
1668           absy = (y < 0) ? -y : y,
1669           kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2,
1670           size, func;
1671
1672    if (C_fitsinbignumhalfdigitp(absy)) {
1673      size = C_fix(C_bignum_size(x));
1674      func = (C_word)bignum_destructive_divide_unsigned_halfdigit;
1675    } else {
1676      size = C_fix(C_bignum_size(x)) + 1; /* Due to normalization */
1677      func = (C_word)bignum_destructive_divide_unsigned_digit;
1678    }
1679
1680    k2 = C_closure(&ka, 7, func, k, x, C_fix(absy),
1681                   /* Return quotient *and* remainder */
1682                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
1683    C_allocate_bignum(3, (C_word)NULL, k2, size, q_negp, r_negp);
1684  }
1685}
1686
1687static void
1688bignum_divrem_fixnum_2(C_word c, C_word self, C_word negated_big)
1689{
1690   C_word k = C_block_item(self, 1);
1691   C_values(4, C_SCHEME_UNDEFINED, k, negated_big, C_fix(0));
1692}
1693
1694void C_ccall
1695C_u_bignum_divrem_bignum(C_word c, C_word self, C_word k, C_word x, C_word y)
1696{
1697  C_word kab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *ka = kab, k2, size,
1698         q_negp = C_mk_bool(C_bignum_negativep(y) ?
1699                            !C_bignum_negativep(x) :
1700                            C_bignum_negativep(x)),
1701         r_negp = C_mk_bool(C_bignum_negativep(x));
1702
1703  switch(bignum_cmp_unsigned(x, y)) {
1704  case 0:
1705    C_values(4, C_SCHEME_UNDEFINED, k,
1706             C_truep(q_negp) ? C_fix(-1) : C_fix(1), C_fix(0));
1707  case -1:
1708    C_values(4, C_SCHEME_UNDEFINED, k, C_fix(0), x);
1709  case 1:
1710    k2 = C_closure(&ka, 9, (C_word)bignum_divide_2_unsigned, k, x, y,
1711                   /* Return quotient *and* remainder */
1712                   C_SCHEME_TRUE, C_SCHEME_TRUE, r_negp,
1713                   /* Will be filled in later */
1714                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
1715    size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
1716    C_allocate_bignum(3, (C_word)NULL, k2, size, q_negp, C_SCHEME_FALSE);
1717  }
1718}
1719
1720void C_ccall
1721C_u_bignum_remainder_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y)
1722{
1723  C_word negp = C_mk_bool(C_bignum_negativep(x));
1724   
1725  y = C_unfix(y);
1726
1727  if (y == 1 || y == -1) {
1728    C_kontinue(k, C_fix(0));
1729  } else if (y == C_MOST_NEGATIVE_FIXNUM) {
1730    /* This is the only case we need to go allocate a bignum for */
1731    C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2, size;
1732
1733    y = C_a_u_i_fix_to_big(&a, C_fix(C_MOST_NEGATIVE_FIXNUM));
1734
1735    /* We can skip bignum_divide_2_unsigned because we need no quotient */
1736    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned_2, k, x, y,
1737                   /* Do not return quotient, do return remainder */
1738                   C_SCHEME_FALSE, C_SCHEME_TRUE, negp,
1739                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
1740    size = C_fix(C_bignum_size(x) + 1);
1741    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1742  } else {
1743    C_word absy = (y < 0) ? -y : y;
1744     
1745    if (C_fitsinbignumhalfdigitp(absy)) {
1746      bignum_destructive_remainder_unsigned_halfdigit(k, x, absy, negp);
1747    } else {
1748      C_word kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2,
1749             size = C_fix(C_bignum_size(x)) + 1; /* Due to normalization */
1750     
1751      k2 = C_closure(&ka, 7, (C_word)bignum_destructive_divide_unsigned_digit,
1752                     k, x, C_fix(absy),
1753                     /* Do not return quotient, do return remainder */
1754                     C_SCHEME_FALSE, C_SCHEME_TRUE, negp);
1755      C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1756    }
1757  }
1758}
1759
1760void C_ccall
1761C_u_bignum_remainder_bignum(C_word c, C_word self, C_word k, C_word x, C_word y)
1762{
1763  C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2, size, negp;
1764
1765  switch(bignum_cmp_unsigned(x, y)) {
1766  case 0:
1767    C_kontinue(k, C_fix(0));
1768  case -1:
1769    C_kontinue(k, x);
1770  case 1:
1771    negp = C_mk_bool(C_bignum_negativep(x));
1772
1773    /* We can skip bignum_divide_2_unsigned because we need no quotient */
1774    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned_2, k, x, y,
1775                   /* Do not return quotient, do return remainder */
1776                   C_SCHEME_FALSE, C_SCHEME_TRUE, negp,
1777                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
1778    size = C_fix(C_bignum_size(x) + 1); /* May need to be normalized */
1779    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1780  }
1781}
1782
1783static void
1784bignum_destructive_remainder_unsigned_halfdigit(C_word k, C_word n, C_word d, C_word negp)
1785{
1786  C_word two_digits,
1787         *start = C_bignum_digits(n),
1788         *scan = start + C_bignum_size(n),
1789         r = 0;
1790
1791  assert((d > 1) && (C_fitsinbignumhalfdigitp(d)));
1792  while (start < scan) {
1793    two_digits = (*--scan);
1794    r = C_BIGNUM_DIGIT_COMBINE(r, C_BIGNUM_DIGIT_HI_HALF(two_digits)) % d;
1795    r = C_BIGNUM_DIGIT_COMBINE(r, C_BIGNUM_DIGIT_LO_HALF(two_digits)) % d;
1796  }
1797  C_kontinue(k, C_truep(negp) ? C_fix(-r) : C_fix(r));
1798}
1799
1800void C_ccall
1801C_basic_quotient(C_word c, C_word self, C_word k, C_word loc, C_word x, C_word y)
1802{
1803  C_word ab[nmax(C_SIZEOF_FLONUM * 2, C_SIZEOF_CLOSURE(5))], *a = ab, k2;
1804
1805  if (x & C_FIXNUM_BIT) {
1806    if (y & C_FIXNUM_BIT) {
1807      C_kontinue(k, C_u_i_fixnum_quotient_checked_loc(loc, x, y));
1808    } else if (C_immediatep(y)) {
1809      C_kontinue(k, C_SCHEME_FALSE);
1810    } else if (C_block_header(y) == C_FLONUM_TAG) {
1811      if (C_isnan(C_flonum_magnitude(y)) ||
1812          C_isinf(C_flonum_magnitude(y)) ||
1813          !C_truep(C_u_i_fpintegerp(y))) {
1814        C_kontinue(k, C_SCHEME_FALSE);
1815      } else {
1816        x = C_a_i_fix_to_flo(&a, 1, x);
1817        C_kontinue(k, C_a_i_flonum_quotient_checked_loc(&a, 3, loc, x, y));
1818      }
1819    } else if (C_IS_BIGNUM_TYPE(y)) {
1820      C_u_integer_quotient(3, (C_word)NULL, k, loc, x, y);
1821    } else {
1822      C_kontinue(k, C_SCHEME_FALSE);
1823    }
1824  } else if (C_immediatep(x)) {
1825    C_kontinue(k, C_SCHEME_FALSE);
1826  } else if (C_block_header(x) == C_FLONUM_TAG) {
1827    if (C_isnan(C_flonum_magnitude(x)) ||
1828        C_isinf(C_flonum_magnitude(x)) ||
1829        !C_truep(C_u_i_fpintegerp(x))) {
1830        C_kontinue(k, C_SCHEME_FALSE);
1831    } else if (y & C_FIXNUM_BIT) {
1832      y = C_a_i_fix_to_flo(&a, 1, y);
1833      C_kontinue(k, C_a_i_flonum_quotient_checked_loc(&a, 3, loc, x, y));
1834    } else if (C_immediatep(y)) {
1835      C_kontinue(k, C_SCHEME_FALSE);
1836    } else if (C_block_header(y) == C_FLONUM_TAG) {
1837      if (C_isnan(C_flonum_magnitude(y)) ||
1838          C_isinf(C_flonum_magnitude(y)) ||
1839          !C_truep(C_u_i_fpintegerp(y))) {
1840        C_kontinue(k, C_SCHEME_FALSE);
1841      } else {
1842        C_kontinue(k, C_a_i_flonum_quotient_checked_loc(&a, 3, loc, x, y));
1843      }
1844    } else if (C_IS_BIGNUM_TYPE(y)) {
1845      if (C_truep(C_u_i_fpintegerp(x))) {
1846        k2 = C_closure(&a, 5, (C_word)quotient_intflo_2,
1847                       k, loc, C_SCHEME_TRUE, y);
1848        C_u_flo_to_int(1, (C_word)NULL, k2, x);
1849      } else { /* May overflow... */
1850        y = C_a_u_i_big_to_flo(&a, 1, y); /* Can't be zero, no need to check */
1851        C_kontinue(k, C_a_i_flonum_quotient(&a, 2, x, y));
1852      }
1853    } else {
1854      C_kontinue(k, C_SCHEME_FALSE);
1855    }
1856  } else if (C_IS_BIGNUM_TYPE(x)) {
1857    if (y & C_FIXNUM_BIT) {
1858      C_u_integer_quotient(3, (C_word)NULL, k, loc, x, y);
1859    } else if (C_immediatep(y)) {
1860      C_kontinue(k, C_SCHEME_FALSE);
1861    } else if (C_block_header(y) == C_FLONUM_TAG) {
1862      if (C_truep(C_u_i_fpintegerp(y))) {
1863        k2 = C_closure(&a, 5, (C_word)quotient_intflo_2,
1864                       k, loc, C_SCHEME_FALSE, x);
1865        C_u_flo_to_int(1, (C_word)NULL, k2, y);
1866      } else if (C_isnan(C_flonum_magnitude(y)) ||
1867                 C_isinf(C_flonum_magnitude(y))) {
1868        C_kontinue(k, C_SCHEME_FALSE);
1869      } else { /* May overflow... */
1870        x = C_a_u_i_big_to_flo(&a, 1, x); /* Can't be zero, no need to check */
1871        C_kontinue(k, C_a_i_flonum_quotient(&a, 2, x, y));
1872      }
1873    } else if (C_IS_BIGNUM_TYPE(y)) {
1874      C_u_bignum_quotient(2, (C_word)NULL, k, x, y);
1875    } else {
1876      C_kontinue(k, C_SCHEME_FALSE);
1877    }
1878  } else {
1879    C_kontinue(k, C_SCHEME_FALSE);
1880  }
1881}
1882
1883static void C_ccall quotient_intflo_2(C_word c, C_word self, C_word intnum)
1884{
1885  C_word k = C_block_item(self, 1),
1886         loc = C_block_item(self, 2),
1887         intnum_is_x =  C_block_item(self, 3),
1888         other_arg = C_block_item(self, 4),
1889         kab[C_SIZEOF_CLOSURE(2)], *ka = kab, k2;
1890 
1891  k2 = C_closure(&ka, 2, (C_word)quotient_intflo_3, k);
1892  if (C_truep(intnum_is_x))
1893    C_u_integer_quotient(3, (C_word)NULL, k2, loc, intnum, other_arg);
1894  else
1895    C_u_integer_quotient(3, (C_word)NULL, k2, loc, other_arg, intnum);
1896}
1897
1898static void C_ccall quotient_intflo_3(C_word c, C_word self, C_word x)
1899{
1900   C_word k = C_block_item(self, 1),
1901          ab[C_SIZEOF_FLONUM], *a = ab;
1902   if (x & C_FIXNUM_BIT) C_kontinue(k, C_a_i_fix_to_flo(&a, 1, x));
1903   else C_kontinue(k, C_a_u_i_big_to_flo(&a, 1, x));
1904}
1905
1906void C_ccall
1907C_u_integer_quotient(C_word c, C_word self, C_word k, C_word loc, C_word x, C_word y)
1908{
1909  if (!(y & C_FIXNUM_BIT)) { /* y is bignum. */
1910    if (x & C_FIXNUM_BIT) {
1911      /* abs(x) < abs(y), so it will always be zero except for this case: */
1912      if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) &&
1913          C_bignum_negated_fitsinfixnump(y)) {
1914        C_kontinue(k, C_fix(-1));
1915      } else {
1916        C_kontinue(k, C_fix(0));
1917      }
1918    } else {
1919      C_u_bignum_quotient(2, (C_word)NULL, k, x, y);
1920    }
1921  } else if (x & C_FIXNUM_BIT) { /* both x and y are fixnum. */
1922    C_kontinue(k, C_u_i_fixnum_quotient_checked_loc(loc, x, y));
1923  } else { /* x is bignum, y is fixnum. */
1924    if (y == C_fix(1)) {
1925      C_kontinue(k, x);
1926    } else if (y == C_fix(0)) {
1927      char *a;
1928      C_div_by_zero_error(C_strloc(a, loc));
1929    } else if (y == C_fix(-1)) {
1930      C_u_bignum_negate(1, (C_word)NULL, k, x);
1931    } else if (y == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
1932      /* This is the only case we need to go allocate a bignum for */
1933      C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2,
1934             negp = C_mk_bool(!C_bignum_negativep(x)), size;
1935
1936      y = C_a_u_i_fix_to_big(&a, C_fix(C_MOST_NEGATIVE_FIXNUM));
1937
1938      k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned, k, x, y,
1939                     /* Return quotient, not remainder */
1940                     C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE,
1941                     /* Will be filled in later */
1942                     C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
1943      size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
1944      C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1945    } else {
1946      C_word negp = C_mk_bool((y < 0) ? !C_bignum_negativep(x) :
1947                                        C_bignum_negativep(x)),
1948             absy = (y < 0) ? -C_unfix(y) : C_unfix(y),
1949             kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2,
1950             size, func;
1951
1952      if (C_fitsinbignumhalfdigitp(absy)) {
1953        size = C_fix(C_bignum_size(x));
1954        func = (C_word)bignum_destructive_divide_unsigned_halfdigit;
1955      } else {
1956        size = C_fix(C_bignum_size(x)) + 1; /* Due to normalization */
1957        func = (C_word)bignum_destructive_divide_unsigned_digit;
1958      }
1959
1960      k2 = C_closure(&ka, 7, func, k, x, C_fix(absy),
1961                     /* Return quotient, not remainder */
1962                     C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
1963      C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1964    }
1965  }
1966}
1967
1968void C_ccall
1969C_u_bignum_quotient(C_word c, C_word self, C_word k, C_word x, C_word y)
1970{
1971  C_word ab[C_SIZEOF_CLOSURE(9)], *a = ab, k2, size,
1972         negp = C_mk_bool(C_bignum_negativep(x) ?
1973                         !C_bignum_negativep(y) :
1974                         C_bignum_negativep(y));
1975
1976  switch(bignum_cmp_unsigned(x, y)) {
1977  case 0:
1978    C_kontinue(k, C_truep(negp) ? C_fix(-1) : C_fix(1));
1979  case -1:
1980    C_kontinue(k, C_fix(0));
1981  case 1:
1982    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned, k, x, y,
1983                   /* Return quotient, not remainder */
1984                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE,
1985                   /* Will be filled in later */
1986                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
1987    size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
1988    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1989  }
1990}
1991
1992static void
1993bignum_destructive_divide_unsigned_halfdigit(C_word c, C_word self, C_word quotient)
1994{
1995  C_word k = C_block_item(self, 1),
1996         numerator = C_block_item(self, 2),
1997         denominator = C_unfix(C_block_item(self, 3)),
1998         /* return_quotient = C_block_item(self, 4), */
1999         return_remainder = C_block_item(self, 5),
2000         remainder_negp = C_block_item(self, 6),
2001         *start = C_bignum_digits(quotient),
2002         *end = start + C_bignum_size(quotient) - 1,
2003         remainder;
2004
2005  bignum_digits_destructive_copy(quotient, numerator);
2006
2007  remainder = bignum_digits_destructive_scale_down(start, end, denominator);
2008  assert(C_fitsinbignumhalfdigitp(remainder));
2009
2010  quotient = C_bignum_normalize(quotient);
2011
2012  if (C_truep(return_remainder)) {
2013    remainder = C_truep(remainder_negp) ? -remainder : remainder;
2014    C_values(4, C_SCHEME_UNDEFINED, k, quotient, C_fix(remainder));
2015  } else {
2016    C_kontinue(k, quotient);
2017  }
2018}
2019
2020/* XXX TODO: This seems unnecessary */
2021static void
2022bignum_destructive_divide_unsigned_digit(C_word c, C_word self, C_word quotient)
2023{
2024  C_word k = C_block_item(self, 1),
2025         numerator = C_block_item(self, 2),
2026         denominator = C_unfix(C_block_item(self, 3)),
2027         return_quotient = C_block_item(self, 4),
2028         return_remainder = C_block_item(self, 5),
2029         remainder_negp = C_block_item(self, 6),
2030         *start = C_bignum_digits(quotient),
2031         *scan = start + C_bignum_size(quotient),
2032         qj, remainder = 0, shift = 0;
2033
2034  /* Because `bignum_divide_digit' requires a normalized denominator. */
2035  while (denominator < ((C_word)1 << (C_BIGNUM_DIGIT_LENGTH-1))) {
2036    denominator <<= 1;
2037    shift++;
2038  }
2039
2040  if (shift) {
2041    bignum_destructive_normalize(quotient, numerator, shift);
2042  } else {
2043    /* Could alo use destructive_normalize here, but this is faster */
2044    bignum_digits_destructive_copy(quotient, numerator);
2045    *(scan-1) = 0; /* Most significant digit, not copied */
2046  }
2047
2048  if (C_truep(return_quotient)) {
2049    while (start < scan) {
2050      remainder = bignum_divide_digit(remainder, (*--scan), denominator, &qj);
2051      (*scan) = qj;
2052    }
2053    quotient = C_bignum_normalize(quotient);
2054  } else {
2055    while (start < scan)
2056      remainder = bignum_divide_digit(remainder, (*--scan), denominator, &qj);
2057  }
2058
2059  if (C_truep(return_remainder)) {
2060    remainder >>= shift;
2061    remainder = C_truep(remainder_negp) ? -remainder : remainder;
2062
2063    if (C_truep(return_quotient))
2064      C_values(4, C_SCHEME_UNDEFINED, k, quotient, C_fix(remainder));
2065    else
2066      C_kontinue(k, C_fix(remainder));
2067  } else {
2068    assert(C_truep(return_quotient));
2069    C_kontinue(k, quotient);
2070  }
2071}
2072
2073static void
2074bignum_destructive_normalize(C_word target, C_word source, C_word shift_left)
2075{
2076  C_word digit, carry = 0,
2077         *scan_source = C_bignum_digits(source),
2078         *scan_target = C_bignum_digits(target),
2079         *end_source = scan_source + C_bignum_size(source),
2080         *end_target = scan_target + C_bignum_size(target),
2081         shift_right = C_BIGNUM_DIGIT_LENGTH - shift_left,
2082         mask = (1L << shift_right) - 1;
2083
2084  while (scan_source < end_source) {
2085    digit = (*scan_source++);
2086    (*scan_target++) = (((digit & mask) << shift_left) | carry);
2087    carry = (digit >> shift_right);
2088  }
2089  if (scan_target < end_target)
2090    (*scan_target) = carry;
2091  else
2092    assert(carry == 0);
2093}
2094
2095/* TODO: This pointer stuff here is suspicious: it's probably slow */
2096static C_word
2097bignum_divide_and_subtract_digit(C_word v1, C_word v2, C_word guess, C_word *u)
2098{
2099  C_word product, diff, sum, carry;
2100
2101  product = v2 * guess;
2102  diff = u[2] - C_BIGNUM_DIGIT_LO_HALF(product);
2103  if (diff < 0) {
2104    u[2] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
2105    carry = C_BIGNUM_DIGIT_HI_HALF(product) + 1;
2106  } else {
2107    u[2] = diff;
2108    carry = C_BIGNUM_DIGIT_HI_HALF(product);
2109  }
2110 
2111  product = v1 * guess + carry;
2112  diff = u[1] - C_BIGNUM_DIGIT_LO_HALF(product);
2113  if (diff < 0) {
2114    u[1] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
2115    carry = C_BIGNUM_DIGIT_HI_HALF(product) + 1;
2116  } else {
2117    u[1] = diff;
2118    carry = C_BIGNUM_DIGIT_HI_HALF(product);
2119    if (carry == 0) return guess; /* DONE */
2120  }
2121
2122  diff = u[0] - carry;
2123  if (diff < 0) {
2124    u[0] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
2125  } else {
2126    u[0] = diff;
2127    return guess; /* DONE */
2128  }
2129
2130  sum = v2 + u[2];
2131  u[2] = sum & C_BIGNUM_HALF_DIGIT_MASK;
2132  carry = C_fitsinbignumhalfdigitp(sum);
2133
2134  sum = v1 + u[1] + carry;
2135  u[1] = sum & C_BIGNUM_HALF_DIGIT_MASK;
2136  carry = C_fitsinbignumhalfdigitp(sum);
2137
2138  u[0] += carry;
2139
2140  return guess - 1;
2141}
2142
2143/* This is a reduced version of the division algorithm, applied to the
2144   case of dividing two bignum digits by one bignum digit.  It is
2145   assumed that the numerator, denominator are normalized. */
2146
2147#define BDD_STEP(qn, j)                                                 \
2148{                                                                       \
2149  uj = u[j];                                                            \
2150  if (uj != v1) {                                                       \
2151    uj_uj1 = C_BIGNUM_DIGIT_COMBINE(uj, u[j + 1]);                      \
2152    guess = uj_uj1 / v1;                                                \
2153    comparand = C_BIGNUM_DIGIT_COMBINE(uj_uj1 % v1, u[j + 2]);          \
2154  } else {                                                              \
2155    guess = C_BIGNUM_HALF_DIGIT_MASK;                                   \
2156    comparand = C_BIGNUM_DIGIT_COMBINE(u[j + 1] + v1, u[j + 2]);        \
2157  }                                                                     \
2158  while ((guess * v2) > comparand) {                                    \
2159    guess -= 1;                                                         \
2160    comparand += v1 << C_BIGNUM_HALF_DIGIT_LENGTH;                      \
2161    if (!C_fitsinbignumdigitp(comparand))                               \
2162      break;                                                            \
2163  }                                                                     \
2164  qn = bignum_divide_and_subtract_digit(v1, v2, guess, &u[j]);          \
2165}
2166
2167/* Because the algorithm from Knuth requires combining the two highest
2168 * digits of u (which we can't fit in a C_word), we instead do two
2169 * such steps, a halfdigit at a time.
2170 */
2171static C_word
2172bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q)
2173{
2174  C_word guess, comparand, v1, v2, uj, uj_uj1, q1, q2, u[4];
2175
2176  if (uh == 0) {
2177    if (ul < v) {
2178      *q = 0;
2179      return ul;
2180    } else if (ul == v) {
2181      *q = 1;
2182      return 0;
2183    }
2184  }
2185  u[0] = C_BIGNUM_DIGIT_HI_HALF(uh);
2186  u[1] = C_BIGNUM_DIGIT_LO_HALF(uh);
2187  u[2] = C_BIGNUM_DIGIT_HI_HALF(ul);
2188  u[3] = C_BIGNUM_DIGIT_LO_HALF(ul);
2189  v1 = C_BIGNUM_DIGIT_HI_HALF(v);
2190  v2 = C_BIGNUM_DIGIT_LO_HALF(v);
2191  BDD_STEP(q1, 0);
2192  BDD_STEP(q2, 1);
2193  *q = C_BIGNUM_DIGIT_COMBINE(q1, q2);
2194  return C_BIGNUM_DIGIT_COMBINE(u[2], u[3]);
2195}
2196
2197#undef BDD_STEP
2198
2199/* Full bignum division */
2200
2201static void
2202bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient)
2203{
2204  C_word numerator = C_block_item(self, 2),
2205         remainder_negp = C_block_item(self, 6),
2206         size = C_fix(C_bignum_size(numerator) + 1);
2207
2208  /* Nice: We can recycle the current closure */
2209  C_set_block_item(self, 0, (C_word)bignum_divide_2_unsigned_2);
2210  C_set_block_item(self, 7, quotient);
2211  C_allocate_bignum(3, (C_word)NULL, self, size, remainder_negp, C_SCHEME_FALSE);
2212}
2213
2214/* For help understanding this algorithm, see:
2215   Knuth, Donald E., "The Art of Computer Programming",
2216   volume 2, "Seminumerical Algorithms"
2217   section 4.3.1, "Multiple-Precision Arithmetic". */
2218
2219static void
2220bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder)
2221{
2222  C_word k = C_block_item(self, 1),
2223         numerator = C_block_item(self, 2),
2224         denominator = C_block_item(self, 3),
2225         return_quotient = C_block_item(self, 4),
2226         return_remainder = C_block_item(self, 5),
2227         /* This one may be overwritten with the remainder */
2228         /* remainder_negp = C_block_item(self, 6), */
2229         quotient = C_block_item(self, 7),
2230         length_d = C_bignum_size(denominator),
2231         d1 = *(C_bignum_digits(denominator) + length_d - 1),
2232         shift = 0;
2233
2234  assert(length_d > 1);
2235  while (d1 < ((C_word)1 << (C_BIGNUM_DIGIT_LENGTH-1))) {
2236    d1 <<= 1;
2237    shift++;
2238  }
2239
2240  if (shift == 0) { /* Already normalized */
2241    bignum_digits_destructive_copy(remainder, numerator);
2242    /* Set most significant digit */
2243    *(C_bignum_digits(remainder) + C_bignum_size(numerator)) = 0;
2244
2245    bignum_destructive_divide_normalized(remainder, denominator, quotient);
2246
2247    if (C_truep(C_and(return_quotient, return_remainder))) {
2248      C_values(4, C_SCHEME_UNDEFINED, k,
2249               C_bignum_normalize(quotient), C_bignum_normalize(remainder));
2250    } else if (C_truep(return_remainder)) {
2251      C_kontinue(k, C_bignum_normalize(remainder));
2252    } else {
2253      assert(C_truep(return_quotient));
2254      C_kontinue(k, C_bignum_normalize(quotient));
2255    }
2256  } else {
2257    /* Requires normalisation of denominator;  Allocate temp bignum for it. */
2258    C_set_block_item(self, 0, (C_word)bignum_divide_2_unsigned_3);
2259    C_set_block_item(self, 6, remainder);
2260    C_set_block_item(self, 8, C_fix(shift));
2261    C_allocate_bignum(3, (C_word)NULL, self, C_fix(length_d),
2262                      C_SCHEME_FALSE, C_SCHEME_FALSE);
2263  }
2264}
2265
2266static void
2267bignum_divide_2_unsigned_3(C_word c, C_word self, C_word tmp_big)
2268{
2269  C_word k = C_block_item(self, 1),
2270         numerator = C_block_item(self, 2),
2271         denominator = C_block_item(self, 3),
2272         return_quotient = C_block_item(self, 4),
2273         return_remainder = C_block_item(self, 5),
2274         remainder = C_block_item(self, 6),
2275         quotient = C_block_item(self, 7),
2276         shift = C_unfix(C_block_item(self, 8));
2277
2278  bignum_destructive_normalize(remainder, numerator, shift);
2279  bignum_destructive_normalize(tmp_big, denominator, shift);
2280  bignum_destructive_divide_normalized(remainder, tmp_big, quotient);
2281
2282  if (C_truep(return_remainder)) {
2283    if (C_truep(return_quotient)) {
2284      C_values(4, C_SCHEME_UNDEFINED, k,
2285               C_bignum_normalize(quotient),
2286               bignum_normalize_shifted(remainder, shift));
2287    } else {
2288      C_kontinue(k, bignum_normalize_shifted(remainder, shift));
2289    }
2290  } else {
2291    assert(C_truep(return_quotient));
2292    C_kontinue(k, C_bignum_normalize(quotient));
2293  }
2294}
2295
2296static void
2297bignum_destructive_divide_normalized(C_word u, C_word v, C_word q)
2298{
2299  C_word u_length = C_bignum_size(u),
2300         v_length = C_bignum_size(v),
2301         *u_start = C_bignum_digits(u),
2302         *u_scan = u_start + u_length,
2303         *u_scan_limit = u_start + v_length,
2304         *u_scan_start = u_scan - v_length,
2305         *v_start = C_bignum_digits(v),
2306         *v_end = v_start + v_length,
2307         *q_scan = (q == C_SCHEME_UNDEFINED) ? NULL :
2308                   (C_bignum_digits(q) + C_bignum_size(q)),
2309         v1 = v_end[-1],
2310         v2 = v_end[-2],
2311         ph, /* high half of double-digit product */
2312         pl, /* low half of double-digit product */
2313         guess, uj, qj,
2314         gh, /* high half-digit of guess */
2315         ch, /* high half of double-digit comparand */
2316         v2l = C_BIGNUM_DIGIT_LO_HALF(v2),
2317         v2h = C_BIGNUM_DIGIT_HI_HALF(v2),
2318         cl, /* low half of double-digit comparand */
2319         gl, /* low half-digit of guess */
2320         gm; /* memory loc for reference parameter */
2321
2322  while (u_scan_limit < u_scan) {
2323    uj = (*--u_scan);
2324    if (uj != v1) {
2325      /* comparand = (((((uj * B) + uj1) % v1) * b) + uj2);
2326         guess = (((uj * B) + uj1) / v1); */
2327      cl = u_scan[-2];
2328      ch = bignum_divide_digit(uj, (u_scan[-1]), v1, (&gm));
2329      guess = gm;
2330    } else {
2331      cl = u_scan[-2];
2332      ch = u_scan[-1] + v1;
2333      guess = C_BIGNUM_DIGIT_MASK;
2334    }
2335    while (1) {
2336      /* product = (guess * v2); */
2337      gl = C_BIGNUM_DIGIT_LO_HALF(guess);
2338      gh = C_BIGNUM_DIGIT_HI_HALF(guess);
2339      pl = v2l * gl;
2340      ph = v2l * gh + v2h * gl + C_BIGNUM_DIGIT_HI_HALF(pl);
2341      pl = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(ph),
2342                                  C_BIGNUM_DIGIT_LO_HALF(pl));
2343      ph = v2h * gh + C_BIGNUM_DIGIT_HI_HALF(ph);
2344      /* if (comparand >= product) */
2345      if ((ch > ph) || ((ch == ph) && (cl >= pl)))
2346        break;
2347      guess--;
2348      /* comparand += (v1 << C_BIGNUM_DIGIT_LENGTH) */
2349      ch += v1;
2350      /* if (comparand >= (B^2)) */
2351      if (!C_fitsinbignumdigitp(ch))
2352        break;
2353    }
2354    qj = bignum_divide_and_subtract(v_start, v_end, guess, (--u_scan_start));
2355    if (q_scan != NULL)
2356      (*--q_scan) = qj;
2357  }
2358}
2359
2360static C_word
2361bignum_divide_and_subtract(C_word *v_start, C_word *v_end, C_word guess, C_word *u_start)
2362{
2363  if (guess == 0) {
2364    return 0;
2365  } else {
2366    C_word *v_scan = v_start,
2367           *u_scan = u_start,
2368           carry = 0,
2369           gl, gh, v, pl, vl, vh, ph, diff, sum;
2370
2371    gl = C_BIGNUM_DIGIT_LO_HALF(guess);
2372    gh = C_BIGNUM_DIGIT_HI_HALF(guess);
2373    while (v_scan < v_end) {
2374      v = (*v_scan++);
2375      vl = C_BIGNUM_DIGIT_LO_HALF(v);
2376      vh = C_BIGNUM_DIGIT_HI_HALF(v);
2377      pl = vl * gl + C_BIGNUM_DIGIT_LO_HALF(carry);
2378      ph = vl * gh + vh * gl + C_BIGNUM_DIGIT_HI_HALF(pl) +
2379                               C_BIGNUM_DIGIT_HI_HALF(carry);
2380      diff = (*u_scan) - C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(ph),
2381                                                C_BIGNUM_DIGIT_LO_HALF(pl));
2382      if (diff < 0) {
2383        (*u_scan++) = diff + ((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
2384        carry = vh * gh + C_BIGNUM_DIGIT_HI_HALF(ph) + 1;
2385      } else {
2386        (*u_scan++) = diff;
2387        carry = vh * gh + C_BIGNUM_DIGIT_HI_HALF(ph);
2388      }
2389    }
2390    if (carry == 0)
2391      return guess;
2392
2393    diff = ((*u_scan) - carry);
2394    if (diff < 0) {
2395      (*u_scan) = diff + ((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
2396    } else {
2397      (*u_scan) = diff;
2398      return guess;
2399    }
2400 
2401    /* Subtraction generated carry, implying guess is one too large.
2402       Add v back in to bring it back down. */
2403    v_scan = v_start;
2404    u_scan = u_start;
2405    carry = 0;
2406    while (v_scan < v_end) {
2407      sum = ((*v_scan++) + (*u_scan) + carry);
2408      if (C_fitsinbignumdigitp(sum)) {
2409        (*u_scan++) = sum;
2410        carry = 0;
2411      } else {
2412        (*u_scan++) = sum & C_BIGNUM_DIGIT_MASK;
2413        carry = 1;
2414      }
2415    }
2416    if (carry) {
2417      sum = (*u_scan) + carry;
2418      (*u_scan) = C_fitsinbignumdigitp(sum) ? sum : (sum  & C_BIGNUM_DIGIT_MASK);
2419    }
2420    return guess - 1;
2421  }
2422}
2423
2424/* Like bignum_normalize, but this also shifts division-normalized numbers */
2425static C_word
2426bignum_normalize_shifted(C_word big, C_word shift_right)
2427{
2428  C_word length = C_bignum_size(big),
2429        *start = C_bignum_digits(big),
2430        *scan = start + length,
2431        digit, carry = 0,
2432        shift_left = C_BIGNUM_DIGIT_LENGTH - shift_right,
2433        mask = (1L << shift_right) - 1;
2434 
2435  while (!(*--scan)) {
2436    if (start == scan) { /* Don't bother with anything else */
2437      return C_fix(0);
2438    }
2439    --length;
2440  }
2441
2442  digit = (*scan);
2443  (*scan) = (digit >> shift_right);
2444  length -= (*scan == 0); /* Add 1 or 0 */
2445  carry = ((digit & mask) << shift_left);
2446 
2447  while (start < scan) {
2448    digit = (*--scan);
2449    (*scan) = ((digit >> shift_right) | carry);
2450    carry = ((digit & mask) << shift_left);
2451  }
2452  assert(carry == 0);
2453  assert(C_bignum_size(big) != length);
2454  assert(length != 1 || *C_bignum_digits(big) != 0);
2455
2456  switch(length) {
2457  case 0:
2458    return C_fix(0);
2459  case 1:
2460    return C_fix(C_bignum_negativep(big) ? -*start : *start);
2461  case 2:
2462    if (C_bignum_negativep(big) && *scan == 1 && *start == 0)
2463      return C_fix(C_MOST_NEGATIVE_FIXNUM);
2464    /* FALLTHROUGH */
2465  default:
2466    /* Mutate vector size of internal bignum vector. */
2467    C_block_header(C_internal_bignum(big)) = (C_STRING_TYPE | C_wordstobytes(length+1));
2468    /* Set internal header. */
2469    C_bignum_header(big) = (C_bignum_header(big) & C_BIGNUM_HEADER_SIGN_BIT) | length;
2470    return big;
2471  }
2472}
2473
Note: See TracBrowser for help on using the repository browser.