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

Last change on this file since 31411 was 31411, checked in by sjamaan, 7 years ago

numbers: Convert quotient fully to core naming conventions

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