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

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

numbers: Convert "remainder" to new style

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