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

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

numbers: Convert remainder procedures to core naming conventions

File size: 97.3 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
1035
1036/**
1037 * Below you will find the functions that have been refactored to
1038 * match the "core" style.
1039 *
1040 * Naming convention idea: Maybe name these fixnum ops differently,
1041 * _p_ for "promoting"?  OTOH, it may be safer and cleaner to rename
1042 * the old fixnum ops to have _fx_ to indicate they run in fixnum mode.
1043 */
1044static void bignum_negate_2(C_word c, C_word self, C_word new_big);
1045static void allocate_bignum_2(C_word c, C_word self, C_word bigvec);
1046static C_word bignum_digits_destructive_scale_up_with_carry(C_word *start, C_word *end, C_word fix_factor, C_word carry);
1047
1048static void bignum_plus_unsigned(C_word k, C_word x, C_word y, C_word negp);
1049static void bignum_plus_unsigned_2(C_word c, C_word self, C_word result);
1050static int bignum_cmp_unsigned(C_word x, C_word y);
1051static void bignum_minus_unsigned(C_word k, C_word x, C_word y);
1052static void bignum_minus_unsigned_2(C_word c, C_word self, C_word result);
1053
1054static void bignum_times_halfdigit_fixnum(C_word k, C_word bigx, C_word fixy, C_word negp);
1055static void bignum_times_halfdigit_fixnum_2(C_word c, C_word self, C_word new_big);
1056static void bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp);
1057static void bignum_times_bignum_unsigned_2(C_word c, C_word self, C_word result);
1058static void digits_to_integer_2(C_word c, C_word self, C_word result);
1059static void bignum_to_digits_2(C_word c, C_word self, C_word working_copy);
1060static void bignum_to_digits_3(C_word c, C_word self, C_word string);
1061static void flo_to_int_2(C_word c, C_word self, C_word result);
1062static C_word ilen(C_uword x);
1063static void bignum_allocate_for_shift(C_word c, C_word self, C_word x);
1064static void bignum_negate_after_shift(C_word c, C_word self, C_word result);
1065static void bignum_actual_shift(C_word c, C_word self, C_word result);
1066static void bignum_abs_2(C_word c, C_word self, C_word result);
1067static void bignum_random_2(C_word c, C_word self, C_word result);
1068static void bignum_maybe_negate_magnitude(C_word k, C_word result);
1069static void bignum_negneg_bitwise_op(C_word c, C_word self, C_word result);
1070static void bignum_posneg_bitwise_op(C_word c, C_word self, C_word result);
1071static void bignum_pospos_bitwise_op(C_word c, C_word self, C_word result);
1072static void bignum_destructive_normalize(C_word target, C_word source, C_word shift_left);
1073static void bignum_destructive_remainder_unsigned_halfdigit(C_word k, C_word n, C_word d, C_word negp);
1074static void bignum_destructive_divide_unsigned_halfdigit(C_word c, C_word self, C_word quotient);
1075static void bignum_destructive_divide_unsigned_digit(C_word c, C_word self, C_word quotient);
1076static C_word bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q);
1077static C_word bignum_divide_and_subtract_digit(C_word v1, C_word v2, C_word guess, C_word *u);
1078static void bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient);
1079static void bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder);
1080static void bignum_divide_2_unsigned_3(C_word c, C_word self, C_word tmp_big);
1081static void bignum_destructive_divide_normalized(C_word u, C_word v, C_word q);
1082static C_word bignum_divide_and_subtract(C_word *v_start, C_word *v_end, C_word guess, C_word *u_start);
1083static C_word bignum_normalize_shifted(C_word bignum, C_word shift_right);
1084
1085/* Copy all the digits from source to target, obliterating what was
1086 * there.  If target is larger than source, the most significant
1087 * digits will remain untouched.
1088 */
1089C_inline void bignum_digits_destructive_copy(C_word target, C_word source)
1090{
1091  C_memcpy(C_bignum_digits(target), C_bignum_digits(source),
1092           /* TODO: This is currently in bytes.  If we change the
1093            * representation that needs to change!
1094            * We subtract the size of the header, too.
1095            */
1096           C_header_size(C_internal_bignum(source))-C_wordstobytes(1));
1097}
1098
1099/* Eventually this will probably need to be integrated into C_2_plus. */
1100void C_ccall
1101C_u_2_fixnum_plus(C_word c, C_word self, C_word k, C_word x, C_word y)
1102{
1103  C_word z, ab[C_SIZEOF_BIGNUM(2)], *a = ab;
1104
1105  /* Exceptional situation: this will cause a real overflow */
1106  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) && y == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
1107    C_kontinue(k, C_bignum2(&a, 1, 0, 2));
1108  } else {
1109    z = C_unfix(x) + C_unfix(y);
1110
1111    /* This code "knows" that both fixnums and bignums have 2 reserved bits */
1112    if(!C_fitsinfixnump(z)) {
1113      /* TODO: function returning either a fixnum or a bignum from a C int */
1114      /* This should help with the C API/FFI too. */
1115      C_kontinue(k, C_bignum2(&a, (z < 0), labs(z) & (C_uword)BIGNUM_DIGIT_MASK, 1));
1116    } else {
1117      C_kontinue(k, C_fix(z));
1118    }
1119  }
1120}
1121
1122void C_ccall
1123C_u_fixnum_plus_bignum(C_word c, C_word self, C_word k, C_word x, C_word y)
1124{
1125  C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab, bigx;
1126  bigx = C_a_u_i_fix_to_big(&a, x);
1127  C_u_2_bignum_plus(2, (C_word)NULL, k, bigx, y);
1128}
1129
1130void C_ccall
1131C_u_2_bignum_plus(C_word c, C_word self, C_word k, C_word x, C_word y)
1132{
1133  if (C_bignum_negativep(x)) {
1134    if (C_bignum_negativep(y)) {
1135      bignum_plus_unsigned(k, x, y, C_SCHEME_TRUE);
1136    } else {
1137      bignum_minus_unsigned(k, y, x);
1138    }
1139  } else {
1140    if (C_bignum_negativep(y)) {
1141      bignum_minus_unsigned(k, x, y);
1142    } else {
1143      bignum_plus_unsigned(k, x, y, C_SCHEME_FALSE);
1144    }
1145  }
1146}
1147
1148void C_ccall
1149C_u_fixnum_minus_bignum(C_word c, C_word self, C_word k, C_word x, C_word y)
1150{
1151  C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab, bigx;
1152  bigx = C_a_u_i_fix_to_big(&a, x);
1153  C_u_2_bignum_minus(2, (C_word)NULL, k, bigx, y);
1154}
1155
1156void C_ccall
1157C_u_bignum_minus_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y)
1158{
1159  C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab, bigy;
1160  bigy = C_a_u_i_fix_to_big(&a, y);
1161  C_u_2_bignum_minus(2, (C_word)NULL, k, x, bigy);
1162}
1163
1164void C_ccall
1165C_u_2_bignum_minus(C_word c, C_word self, C_word k, C_word x, C_word y)
1166{
1167  if (C_bignum_negativep(x)) {
1168    if (C_bignum_negativep(y)) {
1169      bignum_minus_unsigned(k, y, x);
1170    } else {
1171      bignum_plus_unsigned(k, x, y, C_SCHEME_TRUE);
1172    }
1173  } else {
1174    if (C_bignum_negativep(y)) {
1175      bignum_plus_unsigned(k, x, y, C_SCHEME_FALSE);
1176    } else {
1177      bignum_minus_unsigned(k, x, y);
1178    }
1179  }
1180}
1181
1182static void
1183bignum_plus_unsigned(C_word k, C_word x, C_word y, C_word negp)
1184{
1185  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
1186
1187  if (C_bignum_size(y) > C_bignum_size(x)) {
1188    C_word z = x;
1189    x = y;
1190    y = z;
1191  }
1192
1193  k2 = C_closure(&ka, 4, (C_word)bignum_plus_unsigned_2, k, x, y);
1194 
1195  size = C_fix(C_bignum_size(x) + 1); /* One more digit, for possible carry. */
1196  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1197}
1198
1199static void
1200bignum_plus_unsigned_2(C_word c, C_word self, C_word result)
1201{
1202  C_word k = C_block_item(self, 1),
1203         x = C_block_item(self, 2),
1204         y = C_block_item(self, 3),
1205         *scan_x = C_bignum_digits(x),
1206         *end_x = scan_x + C_bignum_size(x),
1207         *scan_y = C_bignum_digits(y),
1208         *end_y = scan_y + C_bignum_size(y),
1209         *scan_r = C_bignum_digits(result),
1210         *end_r = scan_r + C_bignum_size(result),
1211         sum, carry = 0;
1212
1213  /* Move over the two numbers simultaneously, adding digits w/ carry. */
1214  while (scan_y < end_y) {
1215    sum = (*scan_x++) + (*scan_y++) + carry;
1216    if (C_fitsinbignumdigitp(sum)) {
1217      (*scan_r++) = sum;
1218      carry = 0;
1219    } else {
1220      (*scan_r++) = sum & C_BIGNUM_DIGIT_MASK;
1221      carry = 1;
1222    }
1223  }
1224 
1225  /* The end of y, the smaller number.  Propagate carry into the rest of x. */
1226  if (carry) {
1227    while (scan_x < end_x) {
1228      sum = (*scan_x++) + 1;
1229      if (C_fitsinbignumdigitp(sum)) {
1230        (*scan_r++) = sum;
1231        carry = 0;
1232        break;
1233      } else {
1234        (*scan_r++) = sum & C_BIGNUM_DIGIT_MASK;
1235      }
1236    }
1237  }
1238
1239  if (carry) { /* We must've reached the end of x, with still a carry. */
1240    (*scan_r) = 1; /* No need to trim: we're now using the extra digit. */
1241    /* No need to normalize: We started with a bignum and it only grew. */
1242    C_kontinue(k, result);
1243  } else {
1244    /* No more carry.  Copy remaining part of x (if any) to result. */
1245    while (scan_x < end_x)
1246      (*scan_r++) = (*scan_x++);
1247
1248    assert(scan_r == end_r - 1);
1249
1250    *scan_r = 0; /* Ensure trimming works.  TODO: set length directly? */
1251    C_kontinue(k, C_bignum_normalize(result));
1252  }
1253}
1254
1255static int
1256bignum_cmp_unsigned(C_word x, C_word y)
1257{
1258  C_word xlen = C_bignum_size(x), ylen = C_bignum_size(y);
1259
1260  if (xlen < ylen) {
1261    return -1;
1262  } else if (xlen > ylen) {
1263    return 1;
1264  } else if (x == y) {
1265    return 0;
1266  } else {
1267    C_word *startx = C_bignum_digits(x);
1268    C_word *scanx = startx + xlen;
1269    C_word *scany = C_bignum_digits(y) + ylen;
1270
1271    while (startx < scanx) {
1272      C_word xdigit = (*--scanx);
1273      C_word ydigit = (*--scany);
1274      if (xdigit < ydigit)
1275        return -1;
1276      if (xdigit > ydigit)
1277        return 1;
1278    }
1279    return 0;
1280  }
1281}
1282
1283static void
1284bignum_minus_unsigned(C_word k, C_word x, C_word y)
1285{
1286  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size, negp;
1287
1288  switch(bignum_cmp_unsigned(x, y)) {
1289  case 0:             /* x = y, return 0 */
1290    C_kontinue(k, C_fix(0));
1291  case -1:            /* abs(x) < abs(y), return -(abs(y) - abs(x)) */
1292    k2 = C_closure(&ka, 4, (C_word)bignum_minus_unsigned_2, k, y, x);
1293   
1294    size = C_fix(C_bignum_size(y)); /* Maximum size of result is length of y. */
1295    C_allocate_bignum(3, (C_word)NULL, k2, size, C_SCHEME_TRUE, C_SCHEME_FALSE);
1296  case 1:             /* abs(x) > abs(y), return abs(x) - abs(y) */
1297    k2 = C_closure(&ka, 4, (C_word)bignum_minus_unsigned_2, k, x, y);
1298   
1299    size = C_fix(C_bignum_size(x)); /* Maximum size of result is length of x. */
1300    C_allocate_bignum(3, (C_word)NULL, k2, size, C_SCHEME_FALSE, C_SCHEME_FALSE);
1301    break;
1302  }
1303}
1304
1305static void
1306bignum_minus_unsigned_2(C_word c, C_word self, C_word result)
1307{
1308  C_word k = C_block_item(self, 1),
1309         x = C_block_item(self, 2),
1310         y = C_block_item(self, 3),
1311         *scan_x = C_bignum_digits(x),
1312         *end_x = scan_x + C_bignum_size(x),
1313         *scan_r = C_bignum_digits(result),
1314         *scan_y = C_bignum_digits(y),
1315         *end_y = scan_y + C_bignum_size(y),
1316         difference, borrow = 0;
1317
1318  /* Move over the two numbers simultaneously, subtracting digits w/ borrow. */
1319  while (scan_y < end_y) {
1320    difference = (((*scan_x++) - (*scan_y++)) - borrow);
1321    if (difference < 0) {
1322      (*scan_r++) = ((C_word)1 << C_BIGNUM_DIGIT_LENGTH) + difference;
1323      borrow = 1;
1324    } else {
1325      (*scan_r++) = difference;
1326      borrow = 0;
1327    }
1328  }
1329
1330  /* The end of y, the smaller number.  Propagate borrow into the rest of x. */
1331  if (borrow != 0) {
1332    while (scan_x < end_x) {
1333      difference = ((*scan_x++) - borrow);
1334      if (difference < 0) {
1335        /* TODO: Define something like BIGNUM_RADIX if we need it elsewhere */
1336        (*scan_r++) = ((C_word)1 << C_BIGNUM_DIGIT_LENGTH) + difference;
1337      } else {
1338        (*scan_r++) = difference;
1339        borrow = 0;
1340        break;
1341      }
1342    }
1343  }
1344
1345  assert(borrow == 0);
1346
1347  /* Finishing up: Copy remaining part of x (if any) into result. */
1348  while (scan_x < end_x)
1349    (*scan_r++) = (*scan_x++);
1350
1351  C_kontinue(k, C_bignum_normalize(result));
1352}
1353
1354/* TODO: This should probably be renamed C_fixnum_negate to replace
1355 * what's in core.  Unfortunately, that one is allocating inline.
1356 * That one may be renamed to C_u_i_fixnum_negate() or some such.
1357 * TODO: Convert this to be an inline function and move to header?
1358 */
1359void C_ccall
1360C_u_fixnum_neg(C_word c, C_word self, C_word k, C_word x)
1361{
1362  /* Exceptional situation: this will cause an overflow to itself */
1363  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM)) { /* C_fitsinfixnump(x) */
1364    C_word ab[C_SIZEOF_BIGNUM(2)], *a = ab;
1365    C_kontinue(k, C_bignum2(&a, 0, 0, 1));
1366  } else {
1367    C_kontinue(k, C_fix(-C_unfix(x)));
1368  }
1369}
1370
1371C_word C_ccall
1372C_u_i_2_fixnum_gcd(C_word x, C_word y)
1373{
1374   C_word r;
1375   
1376   x = C_unfix(x);
1377   y = C_unfix(y);
1378   
1379   if (x < 0) x = -x;
1380   if (y < 0) y = -y;
1381   
1382   while(y != 0) {
1383     r = x % y;
1384     x = y;
1385     y = r;
1386   }
1387   return C_fix(x);
1388}
1389
1390C_word C_ccall
1391C_a_u_i_2_flonum_gcd(C_word **p, C_word n, C_word x, C_word y)
1392{
1393   double xub, yub, r;
1394   
1395   xub = C_flonum_magnitude(x);
1396   yub = C_flonum_magnitude(y);
1397
1398   if (xub < 0.0) xub = -xub;
1399   if (yub < 0.0) yub = -yub;
1400   
1401   while(yub != 0.0) {
1402     r = fmod(xub, yub);
1403     xub = yub;
1404     yub = r;
1405   }
1406   C_return(C_flonum(p, xub));
1407}
1408
1409void C_ccall
1410C_u_bignum_negate(C_word c, C_word self, C_word k, C_word x)
1411{
1412  C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, negp, size;
1413
1414  negp = C_mk_bool(!C_bignum_negativep(x));
1415  k2 = C_closure(&ka, 3, (C_word)bignum_negate_2, k, x);
1416 
1417  size = C_fix(C_bignum_size(x));
1418  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1419}
1420
1421static void
1422bignum_negate_2(C_word c, C_word self, C_word new_big)
1423{
1424  C_word k = C_block_item(self, 1),
1425         old_big = C_block_item(self, 2);
1426  bignum_digits_destructive_copy(new_big, old_big);
1427  C_kontinue(k, C_bignum_normalize(new_big));
1428}
1429
1430C_word
1431C_u_i_bignum_cmp(C_word x, C_word y)
1432{
1433  if (C_bignum_negativep(x)) {
1434    if (C_bignum_negativep(y)) { /* Largest negative number is smallest */
1435      return C_fix(bignum_cmp_unsigned(y, x));
1436    } else {
1437      return C_fix(-1);
1438    }
1439  } else {
1440    if (C_bignum_negativep(y)) {
1441      return C_fix(1);
1442    } else {
1443      return C_fix(bignum_cmp_unsigned(x, y));
1444    }
1445  }
1446}
1447
1448/* XXX TODO: Maybe pass true/false/bignum as initp, to allow copying data? */
1449void C_ccall
1450C_allocate_bignum(C_word c, C_word self, C_word k, C_word size, C_word negp, C_word initp)
1451{
1452  C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, init;
1453  k2 = C_closure(&ka, 3, (C_word)allocate_bignum_2, k, negp);
1454
1455  init = C_and(initp, C_make_character('\0'));
1456  C_allocate_vector(6, (C_word)NULL, k2,
1457                    C_bytes(C_fixnum_plus(size, C_fix(1))), /* Add header */
1458                    /* Byte vec, initialization, align at 8 bytes */
1459                    C_SCHEME_TRUE, init, C_SCHEME_FALSE);
1460}
1461
1462static void
1463allocate_bignum_2(C_word c, C_word self, C_word bigvec)
1464{
1465  C_word ab[C_SIZEOF_STRUCTURE(2)], *a = ab, bignum,
1466         k = C_block_item(self, 1),
1467         negp = C_truep(C_block_item(self, 2)),
1468         size = C_bytestowords(C_header_size(bigvec))-1;
1469
1470  C_word tagvec = CHICKEN_gc_root_ref(tags);
1471
1472  C_set_block_item(bigvec, 0, negp ? C_BIGNUM_HEADER_SIGN_BIT | size : size);
1473
1474  bignum = C_structure(&a, 2, C_block_item(tagvec, BIG_TAG), bigvec);
1475  C_kontinue(k, bignum);
1476}
1477
1478/* Normalization: scan trailing zeroes, then return a fixnum if the
1479 * value fits, or trim the bignum's length. */
1480C_word C_ccall
1481C_bignum_normalize(C_word big)
1482{
1483  C_word *start = C_bignum_digits(big);
1484  C_word *last_digit = start + C_bignum_size(big) - 1;
1485  C_word *scan = last_digit, length;
1486
1487  while (scan >= start && *scan == 0)
1488    scan--;
1489  length = scan - start + 1;
1490 
1491  switch(length) {
1492  case 0:
1493    return C_fix(0);
1494  case 1:
1495    return C_fix(C_bignum_negativep(big) ? -*start : *start);
1496  case 2:
1497    if (C_bignum_negativep(big) && *scan == 1 && *start == 0)
1498      return C_fix(C_MOST_NEGATIVE_FIXNUM);
1499    /* FALLTHROUGH */
1500  default:
1501    if (scan < last_digit) {
1502      /* Mutate vector size of internal bignum vector. */
1503      C_block_header(C_internal_bignum(big)) = (C_STRING_TYPE | C_wordstobytes(length+1));
1504      /* Set internal header. */
1505      C_bignum_header(big) = (C_bignum_header(big) & C_BIGNUM_HEADER_SIGN_BIT) | length;
1506    }
1507    return big;
1508  }
1509}
1510
1511static C_word
1512bignum_digits_destructive_scale_up_with_carry(C_word *start, C_word *end, C_word factor, C_word carry)
1513{
1514  C_word digit, product_hi, product_lo, *scan = start;
1515
1516  assert(C_fitsinbignumhalfdigitp(carry));
1517
1518  while (scan < end) {
1519    digit = (*scan);
1520    product_lo = factor * C_BIGNUM_DIGIT_LO_HALF(digit) + carry;
1521    product_hi = factor * C_BIGNUM_DIGIT_HI_HALF(digit) +
1522            C_BIGNUM_DIGIT_HI_HALF(product_lo);
1523    (*scan++) = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(product_hi),
1524                                       C_BIGNUM_DIGIT_LO_HALF(product_lo));
1525    carry = C_BIGNUM_DIGIT_HI_HALF(product_hi);
1526  }
1527  return carry;
1528}
1529
1530/* Given (denominator > 1), it is fairly easy to show that
1531   quotient_high fits a bignum halfdigit, after which it is easy to
1532   see that all digits fit a bignum full digit.
1533   
1534   This works because denominator is known to fit a halfdigit, which
1535   means that the remainder modulo denominator will also fit a halfdigit. */
1536static C_word
1537bignum_digits_destructive_scale_down(C_word *start, C_word *end, C_word denominator)
1538{
1539  C_word numerator, remainder = 0, digit, quotient_high, *scan = end;
1540
1541  assert((denominator > 1) && C_fitsinbignumhalfdigitp(denominator));
1542  while (start <= scan) {
1543    digit = *scan;
1544    numerator = C_BIGNUM_DIGIT_COMBINE(remainder, C_BIGNUM_DIGIT_HI_HALF(digit));
1545    quotient_high = (numerator / denominator);
1546    numerator = C_BIGNUM_DIGIT_COMBINE(numerator % denominator,
1547                                       C_BIGNUM_DIGIT_LO_HALF(digit));
1548    (*scan--) = C_BIGNUM_DIGIT_COMBINE(quotient_high, numerator / denominator);
1549    remainder = numerator % denominator;
1550  }
1551  return remainder;
1552}
1553
1554static void
1555bignum_times_halfdigit_fixnum(C_word k, C_word bigx, C_word fixy, C_word negp)
1556{
1557  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
1558
1559  k2 = C_closure(&ka, 4, (C_word)bignum_times_halfdigit_fixnum_2,
1560                 k, bigx, fixy);
1561
1562  size = C_fix(C_bignum_size(bigx) + 1); /* Needs _at most_ one more digit */
1563  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1564}
1565
1566static void
1567bignum_times_halfdigit_fixnum_2(C_word c, C_word self, C_word new_big)
1568{
1569  C_word k = C_block_item(self, 1),
1570         old_bigx = C_block_item(self, 2),
1571         fixy = C_block_item(self, 3),
1572         *digits = C_bignum_digits(new_big),
1573         *end_digit = digits + C_bignum_size(old_bigx);
1574
1575  bignum_digits_destructive_copy(new_big, old_bigx);
1576
1577  /* Scale up, and sanitise the result. TODO: make normalization one op? */
1578  *end_digit = bignum_digits_destructive_scale_up_with_carry(digits, end_digit,
1579                                                             C_unfix(fixy), 0);
1580  C_kontinue(k, C_bignum_normalize(new_big));
1581}
1582
1583void C_ccall
1584C_u_2_fixnum_times(C_word c, C_word self, C_word k, C_word x, C_word y)
1585{
1586  C_word absx, absy, negp;
1587
1588  /* We don't strictly need the abses and negp in all branches... */
1589  absx = C_unfix(x);
1590  absx = absx < 0 ? -absx : absx;
1591  absy = C_unfix(y);
1592  absy = absy < 0 ? -absy : absy;
1593  negp = ((x & C_INT_SIGN_BIT) ? !(y & C_INT_SIGN_BIT) : (y & C_INT_SIGN_BIT));
1594
1595  if (C_fitsinbignumhalfdigitp(absx)) {
1596     if (absx == 0 || absx == 1 || C_fitsinbignumhalfdigitp(absy)) {
1597       C_kontinue(k, C_fix(negp ? -(absx * absy) : (absx * absy)));
1598     } else {
1599       C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab, bigy;
1600       bigy = C_a_u_i_fix_to_big(&a, y);
1601       bignum_times_halfdigit_fixnum(k, bigy, C_fix(absx), C_mk_bool(negp));
1602     }
1603  } else if (C_fitsinbignumhalfdigitp(absy)) {
1604     if (absy == 0 || absy == 1 /*|| C_fitsinbignumhalfdigitp(absx) */) {
1605       C_kontinue(k, C_fix(negp ? -(absx * absy) : (absx * absy)));
1606     } else {
1607       C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab, bigx;
1608       bigx = C_a_u_i_fix_to_big(&a, x);
1609       bignum_times_halfdigit_fixnum(k, bigx, C_fix(absy), C_mk_bool(negp));
1610     }
1611  } else {
1612    C_word ab[C_SIZEOF_FIX_BIGNUM*2], *a = ab, bigx, bigy;
1613    bigx = C_a_u_i_fix_to_big(&a, x);
1614    bigy = C_a_u_i_fix_to_big(&a, y);
1615    bignum_times_bignum_unsigned(k, bigx, bigy, C_mk_bool(negp));
1616  }
1617}
1618
1619void C_ccall
1620C_u_fixnum_times_bignum(C_word c, C_word self, C_word k, C_word x, C_word y)
1621{
1622  C_word negp, absx;
1623 
1624  if (x == C_fix(0))
1625    C_kontinue(k, C_fix(0));
1626  else if (x == C_fix(1))
1627    C_kontinue(k, y);
1628  else if (x == C_fix(-1))
1629    C_u_bignum_negate(1, (C_word)NULL, k, y);
1630
1631  absx = C_unfix(x);
1632  absx = absx < 0 ? -absx : absx;
1633  negp = (x & C_INT_SIGN_BIT) ? !C_bignum_negativep(y) : C_bignum_negativep(y);
1634 
1635  if (C_fitsinbignumhalfdigitp(absx)) {
1636    bignum_times_halfdigit_fixnum(k, y, C_fix(absx), C_mk_bool(negp));
1637  } else {
1638    C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab, bigx, bigy;
1639    bigx = C_a_u_i_fix_to_big(&a, x);
1640    bignum_times_bignum_unsigned(k, bigx, y, C_mk_bool(negp));
1641  }
1642}
1643
1644void C_ccall
1645C_u_2_bignum_times(C_word c, C_word self, C_word k, C_word x, C_word y)
1646{
1647  C_word negp = C_bignum_negativep(x) ?
1648                !C_bignum_negativep(y) :
1649                C_bignum_negativep(y);
1650  bignum_times_bignum_unsigned(k, x, y, C_mk_bool(negp));
1651}
1652
1653/* Multiplication
1654   Maximum value for product_lo or product_hi:
1655        ((R * R) + (R * (R - 2)) + (R - 1))
1656   Maximum value for carry: ((R * (R - 1)) + (R - 1))
1657        where R == 2^HALF_DIGIT_LENGTH */
1658static void
1659bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp)
1660{
1661  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
1662 
1663  if (C_bignum_size(y) > C_bignum_size(x)) { /* Ensure size(x) <= size(y) */
1664    C_word z = x;
1665    x = y;
1666    y = z;
1667  }
1668
1669  k2 = C_closure(&ka, 4, (C_word)bignum_times_bignum_unsigned_2, k, x, y);
1670
1671  size = C_fix(C_bignum_size(x) + C_bignum_size(y));
1672  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
1673}
1674
1675static void
1676bignum_times_bignum_unsigned_2(C_word c, C_word self, C_word result)
1677{
1678  C_word k = C_block_item(self, 1),
1679         x = C_block_item(self, 2),
1680         y = C_block_item(self, 3),
1681
1682         carry, y_digit_lo, y_digit_hi, x_digit_lo,
1683         x_digit_hi, product_lo, *scan_r, *scan_y,
1684         x_digit, y_digit, product_hi,
1685         *scan_x = C_bignum_digits(x),
1686         *end_x = scan_x + C_bignum_size(x),
1687         *start_y = C_bignum_digits(y),
1688         *end_y = start_y + C_bignum_size(y),
1689         *start_r = C_bignum_digits(result);
1690
1691  while (scan_x < end_x) {
1692    x_digit = (*scan_x++);
1693    x_digit_lo = C_BIGNUM_DIGIT_LO_HALF(x_digit);
1694    x_digit_hi = C_BIGNUM_DIGIT_HI_HALF(x_digit);
1695    carry = 0;
1696    scan_y = start_y;
1697    scan_r = (start_r++);
1698
1699    while (scan_y < end_y) {
1700      y_digit = (*scan_y++);
1701      y_digit_lo = C_BIGNUM_DIGIT_LO_HALF(y_digit);
1702      y_digit_hi = C_BIGNUM_DIGIT_HI_HALF(y_digit);
1703
1704      product_lo = (*scan_r) +
1705                   x_digit_lo * y_digit_lo +
1706                   C_BIGNUM_DIGIT_LO_HALF(carry);
1707
1708      product_hi = x_digit_hi * y_digit_lo +
1709                   x_digit_lo * y_digit_hi +
1710                   C_BIGNUM_DIGIT_HI_HALF(product_lo) +
1711                   C_BIGNUM_DIGIT_HI_HALF(carry);
1712
1713      (*scan_r++) = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(product_hi),
1714                                           C_BIGNUM_DIGIT_LO_HALF(product_lo));
1715
1716      carry = x_digit_hi * y_digit_hi + C_BIGNUM_DIGIT_HI_HALF(product_hi);
1717    }
1718    (*scan_r) += carry;
1719  }
1720  C_kontinue(k, C_bignum_normalize(result));
1721}
1722
1723void C_ccall
1724C_digits_to_integer(C_word c, C_word self, C_word k, C_word str,
1725                    C_word start, C_word end, C_word radix, C_word negp)
1726{
1727  C_word kab[C_SIZEOF_CLOSURE(6)], *ka = kab, k2, size;
1728  size_t nbits;
1729
1730  assert((C_unfix(radix) > 1) && C_fitsinbignumhalfdigitp(C_unfix(radix)));
1731 
1732  if (start == end) {
1733    C_kontinue(k, C_SCHEME_FALSE);
1734  } else {
1735    k2 = C_closure(&ka, 6, (C_word)digits_to_integer_2, k, str, start, end, radix);
1736 
1737    nbits = (C_unfix(end) - C_unfix(start)) * ilen(C_unfix(radix));
1738    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(nbits));
1739    /* XXX: Why initialize? */
1740    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
1741  }
1742}
1743
1744/* XXX TODO: This should be faster, dammit! */
1745static void
1746digits_to_integer_2(C_word c, C_word self, C_word result)
1747{
1748  C_word k = C_block_item(self, 1),
1749         str = C_block_item(self, 2),
1750         start = C_unfix(C_block_item(self, 3)),
1751         end = C_unfix(C_block_item(self, 4)),
1752         radix = C_unfix(C_block_item(self, 5)),
1753         first_digit = C_unfix(C_block_item(self, 6)),
1754         *digits = C_bignum_digits(result),
1755         *last_digit = digits, /* Initially, bignum is all zeroes */
1756         big_digit = 0, factor = radix,
1757         next_big_digit, next_factor,
1758         carry, str_digit;
1759  char *str_scan = C_c_string(str) + start,
1760       *str_end = C_c_string(str) + end;
1761
1762  /* Hash characters in numbers are mapped to 0 */
1763# define HEXDIGIT_CHAR_TO_INT(x)                                        \
1764    (((x) == '#') ? 0 :                                                 \
1765     (((x) >= (int)'a') ? ((x) - (int)'a' + 10) : ((x) - (int)'0')))
1766
1767  /* This tries to save up as much as possible in the local C_word
1768   * big_digit, and only when it exceeds what we would be able to
1769   * multiply easily, we scale up the bignum and add what we saved up.
1770   */
1771  while (str_scan < str_end) {
1772    str_digit = HEXDIGIT_CHAR_TO_INT(C_tolower((int)*str_scan));
1773    str_scan++;
1774
1775    next_big_digit = big_digit * radix;
1776    next_big_digit += str_digit;
1777    next_factor = factor * radix;
1778
1779    if (str_digit >= radix || str_digit < 0) {
1780      C_kontinue(k, C_SCHEME_FALSE);
1781    } else if (C_fitsinbignumhalfdigitp(next_big_digit) &&
1782               C_fitsinbignumhalfdigitp(next_factor)) {
1783      factor = next_factor;
1784      big_digit = next_big_digit;
1785    } else {
1786      carry = bignum_digits_destructive_scale_up_with_carry(
1787              digits, last_digit, factor, big_digit);
1788
1789      if (carry) (*last_digit++) = carry; /* Move end */
1790
1791      big_digit = str_digit;
1792      factor = radix;
1793    }
1794  }
1795# undef HEXDIGIT_CHAR_TO_INT
1796
1797  /* Final step.  We always must do this, because the loop never
1798   * processes the "current" character into the bignum (lookahead 1).
1799   */
1800  carry = bignum_digits_destructive_scale_up_with_carry(
1801          digits, last_digit, factor, big_digit);
1802  if (carry) (*last_digit++) = carry; /* Move end */
1803
1804  C_kontinue(k, C_bignum_normalize(result));
1805}
1806
1807void C_ccall
1808C_u_bignum_to_digits(C_word c, C_word self, C_word k, C_word value, C_word radix)
1809{
1810  /* This copies the bignum over into a working copy that can be mutated. */
1811  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size,
1812         negp = C_mk_bool(C_bignum_negativep(value));
1813
1814  k2 = C_closure(&ka, 4, (C_word)bignum_to_digits_2, k, value, radix);
1815
1816  size = C_fix(C_bignum_size(value));
1817  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1818}
1819
1820static void
1821bignum_to_digits_2(C_word c, C_word self, C_word working_copy)
1822{
1823  C_word k = C_block_item(self, 1),
1824         old_big = C_block_item(self, 2),
1825         radix = C_unfix(C_block_item(self, 3)),
1826         len;
1827
1828  bignum_digits_destructive_copy(working_copy, old_big);
1829
1830  assert(radix > 1 && C_fitsinbignumhalfdigitp(radix));
1831
1832  /* A sloppy over-approximation of the number of radix digits we
1833   * need.  Due to the algorithm chopping off slices of "base" size at
1834   * a time, we add one to the bignum size, so we are sure we can
1835   * handle any leading zeroes.  The extra 1 we add at the end is for
1836   * an optional sign.
1837   */
1838  len = (C_bignum_size(old_big)+1)*C_BIGNUM_DIGIT_LENGTH / (ilen(radix)-1) + 1;
1839
1840  /* Nice: We can recycle the current closure */
1841  C_set_block_item(self, 0, (C_word)bignum_to_digits_3);
1842  /* item 1 is still the same continuation */
1843  C_set_block_item(self, 2, working_copy);
1844
1845  C_allocate_vector(6, (C_word)NULL, self, C_fix(len),
1846                    /* Byte vec, no initialization, align at 8 bytes */
1847                    C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
1848}
1849
1850/* XXX TODO: This should be MUCH faster, dammit! */
1851static void
1852bignum_to_digits_3(C_word c, C_word self, C_word string)
1853{
1854  static char *characters = "0123456789abcdef";
1855  C_word k = C_block_item(self, 1),
1856         working_copy = C_block_item(self, 2),
1857         radix = C_unfix(C_block_item(self, 3)),
1858         *start = C_bignum_digits(working_copy),
1859         *scan = start + C_bignum_size(working_copy) - 1,
1860         len = C_header_size(string),
1861         digit, steps, i, base;
1862  char *buf = C_c_string(string), *index = buf + len - 1;
1863
1864  /* Calculate the largest power of radix that fits a halfdigit:
1865   * steps = log10(2^halfdigit_bits), base = 10^steps
1866   */
1867  for(steps = 0, base = radix; C_fitsinbignumhalfdigitp(base); base *= radix)
1868    steps++;
1869
1870  base /= radix; /* Back down: we overshot in the loop */
1871
1872  while (start <= scan) {
1873    digit = bignum_digits_destructive_scale_down(start, scan, base);
1874
1875    if (*scan == 0) scan--; /* Adjust if we exhausted the highest digit */
1876
1877    for(i = 0; i < steps; ++i) {
1878      *index-- = characters[digit % radix];
1879      digit /= radix;
1880    }
1881  }
1882  assert(index >= buf-1);
1883
1884  /* Move index onto first nonzero digit (advance, and skip leading zeroes) */
1885  while(*++index == '0');
1886
1887  if (C_bignum_negativep(working_copy))
1888    *--index = '-';
1889
1890  len -= index - buf; /* Shorten with distance between start and index. */
1891  C_memmove(buf, index, len); /* Move start of number to begin of string. */
1892  C_block_header(string) = C_STRING_TYPE | len; /* Mutate string length. */
1893  C_kontinue(k, string);
1894}
1895
1896C_word
1897C_a_u_i_big_to_flo(C_word **p, C_word n, C_word bignum)
1898{
1899  double accumulator = 0;
1900  C_word *start = C_bignum_digits(bignum);
1901  C_word *scan = start + C_bignum_size(bignum);
1902  while (start < scan) {
1903    accumulator *= (C_word)1 << C_BIGNUM_DIGIT_LENGTH;
1904    accumulator += (*--scan);
1905  }
1906  return C_flonum(p, (C_bignum_negativep(bignum) ? -accumulator : accumulator));
1907}
1908
1909void C_ccall
1910C_u_flo_to_int(C_word c, C_word self, C_word k, C_word x)
1911{
1912  int exponent;
1913  double significand = frexp(C_flonum_magnitude(x), &exponent);
1914
1915  if (exponent <= 0) {
1916    C_kontinue(k, C_fix(0));
1917  } else if (exponent == 1) { /* TODO: check significand * 2^exp fits fixnum? */
1918    C_kontinue(k, significand < 0.0 ? C_fix(-1) : C_fix(1));
1919  } else {
1920    C_word kab[C_SIZEOF_CLOSURE(4) + C_SIZEOF_FLONUM], *ka = kab, k2, size,
1921           negp = C_mk_bool(C_flonum_magnitude(x) < 0.0),
1922           sign = C_flonum(&ka, fabs(significand));
1923
1924    k2 = C_closure(&ka, 4, (C_word)flo_to_int_2, k, C_fix(exponent), sign);
1925
1926    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(exponent));
1927    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1928  }
1929}
1930
1931static void
1932flo_to_int_2(C_word c, C_word self, C_word result)
1933{
1934  C_word digit, k = C_block_item(self, 1),
1935         exponent = C_unfix(C_block_item(self, 2)),
1936         odd_bits = exponent % C_BIGNUM_DIGIT_LENGTH,
1937         *start = C_bignum_digits(result),
1938         *scan = start + C_bignum_size(result);
1939  /* It's always true that 0.5 <= s < 1 */
1940  double significand = C_flonum_magnitude(C_block_item(self, 3));
1941
1942  if (odd_bits > 0) { /* Handle most significant digit first */
1943    significand *= (C_word)1 << odd_bits;
1944    digit = (C_word)significand;
1945    (*--scan) = digit;
1946    significand -= (double)digit;
1947  }
1948
1949  while (start < scan && significand > 0) {
1950    significand *= (C_word)1 << C_BIGNUM_DIGIT_LENGTH;
1951    digit = (C_word)significand;
1952    (*--scan) = digit;
1953    significand -= (double)digit;
1954  }
1955
1956  /* Finish up by clearing any remaining, lower, digits */
1957  while (start < scan)
1958    (*--scan) = 0;
1959
1960  C_kontinue(k, C_bignum_normalize(result));
1961}
1962
1963/*
1964 * From Hacker's Delight by Henry S. Warren
1965 * based on a modified nlz() from section 5-3 (fig. 5-7)
1966 */
1967static C_word
1968ilen(C_uword x)
1969{
1970  C_uword y;
1971  C_word n = 0;
1972
1973#ifdef C_SIXTY_FOUR
1974  y = x >> 32; if (y != 0) { n += 32; x = y; }
1975#endif
1976  y = x >> 16; if (y != 0) { n += 16; x = y; }
1977  y = x >>  8; if (y != 0) { n +=  8; x = y; }
1978  y = x >>  4; if (y != 0) { n +=  4; x = y; }
1979  y = x >>  2; if (y != 0) { n +=  2; x = y; }
1980  y = x >>  1; if (y != 0) return n + 2;
1981  return n + x;
1982}
1983
1984C_word
1985C_u_i_int_length(C_word x)
1986{
1987  if (x & C_FIXNUM_BIT) {
1988    x = C_unfix(x);
1989    return C_fix(ilen((x < 0) ? ~x : x));
1990  } else { /* bignum */
1991    C_word len_1 = C_bignum_size(x) - 1,
1992           result = len_1 * BIGNUM_DIGIT_LENGTH,
1993           *startx = C_bignum_digits(x),
1994           *last_digit = C_bignum_digits(x) + len_1,
1995           last_digit_length = ilen(*last_digit);
1996
1997    /* If *only* the highest bit is set, negating results in one less bit */
1998    if (C_bignum_negativep(x) && *last_digit == (1 << (last_digit_length-1))) {
1999      while(startx < last_digit && *startx == 0) ++startx;
2000      if (startx == last_digit) result--;
2001    }
2002    return C_fix(result + last_digit_length);
2003  }
2004}
2005
2006void C_ccall
2007C_u_int_shift_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
2008{
2009  C_word kab[C_SIZEOF_FIX_BIGNUM + C_SIZEOF_CLOSURE(3) + C_SIZEOF_CLOSURE(2)],
2010         *ka = kab, k2, k3, size;
2011
2012  if (y == C_fix(0) || x == C_fix(0)) { /* Done (no shift) */
2013    C_kontinue(k, x);
2014  } else if (x & C_FIXNUM_BIT) {
2015    /* TODO: This should probably see if shifting fits a fixnum, first */
2016    x = C_a_u_i_fix_to_big(&ka, x);
2017  }
2018
2019  /* Invert all the bits before shifting right a negative value */
2020  if (C_bignum_negativep(x) && C_unfix(y) < 0) {
2021    /* When done shifting, invert again */
2022    k3 = C_closure(&ka, 2, (C_word)bignum_negate_after_shift, k);
2023    /* Before shifting, allocate the bignum */
2024    k2 = C_closure(&ka, 3, (C_word)bignum_allocate_for_shift, k3, y);
2025    /* Actually invert by subtracting: -1 - x */
2026    C_u_fixnum_minus_bignum(2, (C_word)NULL, k2, C_fix(-1), x);
2027  } else {
2028    k2 = C_closure(&ka, 3, (C_word)bignum_allocate_for_shift, k, y);
2029    C_kontinue(k2, x);
2030  }
2031}
2032
2033static void
2034bignum_allocate_for_shift(C_word c, C_word self, C_word x)
2035{
2036  C_word k = C_block_item(self, 1),
2037         y = C_block_item(self, 2),
2038         uy = C_unfix(y),
2039         negp, digit_offset, bit_offset,
2040         ab[C_SIZEOF_FIX_BIGNUM + C_SIZEOF_CLOSURE(6)], *a = ab, k2, size;
2041
2042  if (x & C_FIXNUM_BIT) /* Normalisation may happen after negation */
2043    x = C_a_u_i_fix_to_big(&a, x);
2044
2045  negp = C_mk_bool(C_bignum_negativep(x));
2046 
2047  /* uy is guaranteed not to be 0 here */
2048  if (uy > 0) {
2049    digit_offset = uy / C_BIGNUM_DIGIT_LENGTH;
2050    bit_offset =   uy % C_BIGNUM_DIGIT_LENGTH;
2051
2052    k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k,
2053                   x, C_SCHEME_TRUE, C_fix(digit_offset), C_fix(bit_offset));
2054    size = C_fix(C_bignum_size(x) + digit_offset + 1);
2055    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
2056  } else if (-uy >= C_bignum_size(x) * (C_word)C_BIGNUM_DIGIT_LENGTH) {
2057    /* All bits are shifted out, just return 0 */
2058    C_kontinue(k, C_fix(0));
2059  } else {
2060    digit_offset = -uy / BIGNUM_DIGIT_LENGTH;
2061    bit_offset =   -uy % BIGNUM_DIGIT_LENGTH;
2062   
2063    k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k,
2064                   x, C_SCHEME_FALSE, C_fix(digit_offset), C_fix(bit_offset));
2065    size = C_fix(C_bignum_size(x) - digit_offset);
2066    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
2067  }
2068}
2069
2070static void
2071bignum_negate_after_shift(C_word c, C_word self, C_word result)
2072{
2073  C_word k = C_block_item(self, 1),
2074         ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
2075  if (result & C_FIXNUM_BIT) /* Normalisation may happen after shift */
2076    C_kontinue(k, C_fix(-1 - C_unfix(result)));
2077  else
2078    C_u_fixnum_minus_bignum(2, (C_word)NULL, k, C_fix(-1), result);
2079}
2080
2081static void
2082bignum_actual_shift(C_word c, C_word self, C_word result)
2083{
2084  C_word k = C_block_item(self, 1),
2085         x = C_block_item(self, 2),
2086         shift_left = C_truep(C_block_item(self, 3)),
2087         digit_offset = C_unfix(C_block_item(self, 4)),
2088         bit_offset = C_unfix(C_block_item(self, 5)),
2089         *scanx, *scanr, *end;
2090
2091  if (shift_left) {
2092    scanr = C_bignum_digits(result) + digit_offset;
2093    scanx = C_bignum_digits(x);
2094    end = scanx + C_bignum_size(x);
2095   
2096    while (scanx < end) {
2097      *scanr = *scanr | (*scanx & C_BIGNUM_DIGIT_MASK) << bit_offset;
2098      *scanr = *scanr & C_BIGNUM_DIGIT_MASK;
2099      scanr++;
2100      *scanr = *scanx++ >> (C_BIGNUM_DIGIT_LENGTH - bit_offset);
2101      *scanr = *scanr & C_BIGNUM_DIGIT_MASK;
2102    }
2103  } else {
2104    scanr = C_bignum_digits(result);
2105    scanx = C_bignum_digits(x) + digit_offset;
2106    end = scanr + C_bignum_size(result) - 1;
2107   
2108    while (scanr < end) {
2109      *scanr =  (*scanx++ & C_BIGNUM_DIGIT_MASK) >> bit_offset;
2110      *scanr = (*scanr | 
2111        *scanx << (C_BIGNUM_DIGIT_LENGTH - bit_offset)) & C_BIGNUM_DIGIT_MASK;
2112      scanr++;
2113    }
2114    *scanr =  (*scanx++ & C_BIGNUM_DIGIT_MASK) >> bit_offset;
2115  }
2116  C_kontinue(k, C_bignum_normalize(result));
2117}
2118
2119void C_ccall
2120C_u_bignum_abs(C_word c, C_word self, C_word k, C_word big)
2121{
2122  if (!C_bignum_negativep(big)) {
2123    C_kontinue(k, big);
2124  } else {
2125    C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, size;
2126
2127    k2 = C_closure(&ka, 3, (C_word)bignum_abs_2, k, big);
2128
2129    size = C_fix(C_bignum_size(big));
2130    C_allocate_bignum(3, (C_word)NULL, k2, size, C_SCHEME_FALSE, C_SCHEME_FALSE);
2131  }
2132}
2133
2134static void
2135bignum_abs_2(C_word c, C_word self, C_word result)
2136{
2137  C_word k = C_block_item(self, 1),
2138         old_big = C_block_item(self, 2);
2139
2140  /* Just copy the old bignum's digits as-is */
2141  bignum_digits_destructive_copy(result, old_big);
2142  C_kontinue(k, result);
2143}
2144
2145/*
2146 * This random number generator is very simple. Probably too simple...
2147 */
2148C_word C_ccall
2149C_u_i_bignum_randomize(C_word bignum)
2150{
2151  C_word seed = 0,
2152         *scan = C_bignum_digits(bignum),
2153         *end = scan + C_bignum_size(bignum);
2154
2155  /* What a cheap way to initialize the random generator. I feel dirty! */
2156  while (scan < end)
2157    seed ^= *scan++;
2158
2159  srand(seed);
2160  return C_SCHEME_UNDEFINED;
2161}
2162
2163void C_ccall
2164C_u_bignum_random(C_word c, C_word self, C_word k, C_word max)
2165{
2166  C_word k2, kab[C_SIZEOF_CLOSURE(4)], *ka = kab, size,
2167         max_len, max_bits, max_top_digit, d, negp;
2168
2169  max_len = C_bignum_size(max);
2170  max_top_digit = d = C_bignum_digits(max)[max_len - 1];
2171 
2172  max_bits = (max_len - 1) * C_BIGNUM_DIGIT_LENGTH;
2173  while(d) {
2174    max_bits++;
2175    d >>= 1;
2176  }
2177  /* Subtract/add one because we don't want zero to be over-represented */
2178  size = ((double)rand())/(RAND_MAX + 1.0) * (double)(max_bits - 1);
2179  size = C_fix(C_BIGNUM_BITS_TO_DIGITS(size + 1));
2180
2181  negp = C_mk_bool(C_bignum_negativep(max));
2182  k2 = C_closure(&ka, 4, (C_word)bignum_random_2, k, C_fix(max_top_digit), C_fix(max_len));
2183  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2184}
2185
2186static void
2187bignum_random_2(C_word c, C_word self, C_word result)
2188{
2189  C_word k = C_block_item(self, 1),
2190         max_top_digit = C_unfix(C_block_item(self, 2)),
2191         max_len = C_unfix(C_block_item(self, 3)),
2192         *scan = C_bignum_digits(result),
2193         *end = scan + C_bignum_size(result); /* Go to just before the end. */
2194
2195  while(scan < end)
2196    *scan++ = ((double)rand())/(RAND_MAX + 1.0) * (double)((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
2197  /*
2198   * Last word is special when length is max_len: It must be less than
2199   * max's most significant digit, instead of BIGNUM_RADIX.
2200   */
2201  if (max_len == C_bignum_size(result))
2202    *scan = ((double)rand())/(RAND_MAX + 1.0) * (double)max_top_digit;
2203  else
2204    *scan = ((double)rand())/(RAND_MAX + 1.0) * (double)((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
2205
2206  C_kontinue(k, C_bignum_normalize(result));
2207}
2208
2209
2210/* op identifies the operator: 0 means AND, 1 means IOR, 2 means XOR */
2211void C_ccall
2212C_u_2_integer_bitwise_op(C_word c, C_word self, C_word k, C_word op, C_word x, C_word y)
2213{
2214  C_word kab[C_SIZEOF_FIX_BIGNUM*2 + C_SIZEOF_CLOSURE(5)], *ka = kab, k2,
2215         size, negp;
2216 
2217  if (x & C_FIXNUM_BIT)
2218    x = C_a_u_i_fix_to_big(&ka, x);
2219
2220  if (y & C_FIXNUM_BIT)
2221    y = C_a_u_i_fix_to_big(&ka, y);
2222
2223# define nmax(x, y)     ((x) > (y) ? (x) : (y)) /* From runtime.c */
2224
2225  if (C_bignum_negativep(x)) {
2226    if (C_bignum_negativep(y)) {
2227      negp = C_mk_bool(op == C_fix(0) || op == C_fix(1)); /* and / ior */
2228      size = C_fix(nmax(C_bignum_size(x) + 1, C_bignum_size(y) + 1));
2229      k2 = C_closure(&ka, 5, (C_word)bignum_negneg_bitwise_op, k, op, x, y);
2230    } else {
2231      negp = C_mk_bool(op == C_fix(1) || op == C_fix(2)); /* ior / xor */
2232      size = C_fix(nmax(C_bignum_size(y), C_bignum_size(x) + 1)); /*!*/
2233      k2 = C_closure(&ka, 5, (C_word)bignum_posneg_bitwise_op, k, op, y, x);
2234    }
2235  } else {
2236    if (C_bignum_negativep(y)) {
2237      negp = C_mk_bool(op == C_fix(1) || op == C_fix(2)); /* ior / xor */
2238      size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y) + 1));
2239      k2 = C_closure(&ka, 5, (C_word)bignum_posneg_bitwise_op, k, op, x, y);
2240    } else {
2241      negp = C_SCHEME_FALSE;
2242      size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y)));
2243      k2 = C_closure(&ka, 5, (C_word)bignum_pospos_bitwise_op, k, op, x, y);
2244    }
2245  }
2246
2247# undef nmax
2248
2249  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2250}
2251
2252static void
2253bignum_pospos_bitwise_op(C_word c, C_word self, C_word result)
2254{
2255  C_word k = C_block_item(self, 1),
2256         op = C_block_item(self, 2),
2257         arg1 = C_block_item(self, 3),
2258         arg2 = C_block_item(self, 4),
2259         *scanr = C_bignum_digits(result),
2260         *endr = scanr + C_bignum_size(result),
2261         *scan1 = C_bignum_digits(arg1),
2262         *end1 = scan1 + C_bignum_size(arg1),
2263         *scan2 = C_bignum_digits(arg2),
2264         *end2 = scan2 + C_bignum_size(arg2),
2265         digit1, digit2;
2266
2267  while (scanr < endr) {
2268    digit1 = (scan1 < end1) ? *scan1++ : 0;
2269    digit2 = (scan2 < end2) ? *scan2++ : 0;
2270    *scanr++ = (op == C_fix(0)) ? digit1 & digit2 :
2271               (op == C_fix(1)) ? digit1 | digit2 :
2272                                  digit1 ^ digit2;
2273  }
2274  C_kontinue(k, C_bignum_normalize(result));
2275}
2276
2277static void
2278bignum_posneg_bitwise_op(C_word c, C_word self, C_word result)
2279{
2280  C_word k = C_block_item(self, 1),
2281         op = C_block_item(self, 2),
2282         arg1 = C_block_item(self, 3),
2283         arg2 = C_block_item(self, 4),
2284         *scanr = C_bignum_digits(result),
2285         *endr = scanr + C_bignum_size(result),
2286         *scan1 = C_bignum_digits(arg1),
2287         *end1 = scan1 + C_bignum_size(arg1),
2288         *scan2 = C_bignum_digits(arg2),
2289         *end2 = scan2 + C_bignum_size(arg2),
2290         digit1, digit2, carry2 = 1;
2291
2292  while (scanr < endr) {
2293    digit1 = (scan1 < end1) ? *scan1++ : 0;
2294    digit2 = (~((scan2 < end2) ? *scan2++ : 0) & C_BIGNUM_DIGIT_MASK) + carry2;
2295
2296    if (C_fitsinbignumdigitp(digit2)) {
2297      carry2 = 0;
2298    } else {
2299      digit2 &= C_BIGNUM_DIGIT_MASK;
2300      carry2 = 1;
2301    }
2302   
2303    *scanr++ = (op == C_fix(0)) ? digit1 & digit2 :
2304               (op == C_fix(1)) ? digit1 | digit2 :
2305                                  digit1 ^ digit2;
2306  }
2307
2308  bignum_maybe_negate_magnitude(k, result);
2309}
2310
2311static void
2312bignum_negneg_bitwise_op(C_word c, C_word self, C_word result)
2313{
2314  C_word k = C_block_item(self, 1),
2315         op = C_block_item(self, 2),
2316         arg1 = C_block_item(self, 3),
2317         arg2 = C_block_item(self, 4),
2318         *scanr = C_bignum_digits(result),
2319         *endr = scanr + C_bignum_size(result),
2320         *scan1 = C_bignum_digits(arg1),
2321         *end1 = scan1 + C_bignum_size(arg1),
2322         *scan2 = C_bignum_digits(arg2),
2323         *end2 = scan2 + C_bignum_size(arg2),
2324         digit1, digit2, carry1 = 1, carry2 = 1;
2325
2326  while (scanr < endr) {
2327    digit1 = (~((scan1 < end1) ? *scan1++ : 0) & C_BIGNUM_DIGIT_MASK) + carry1;
2328    digit2 = (~((scan2 < end2) ? *scan2++ : 0) & C_BIGNUM_DIGIT_MASK) + carry2;
2329
2330    if (C_fitsinbignumdigitp(digit1)) {
2331      carry1 = 0;
2332    } else {
2333      digit1 &= C_BIGNUM_DIGIT_MASK;
2334      carry1 = 1;
2335    }
2336    if (C_fitsinbignumdigitp(digit2)) {
2337      carry2 = 0;
2338    } else {
2339      digit2 &= C_BIGNUM_DIGIT_MASK;
2340      carry2 = 1;
2341    }
2342   
2343    *scanr++ = (op == C_fix(0)) ? digit1 & digit2 :
2344               (op == C_fix(1)) ? digit1 | digit2 :
2345                                  digit1 ^ digit2;
2346  }
2347
2348  bignum_maybe_negate_magnitude(k, result);
2349}
2350
2351static void
2352bignum_maybe_negate_magnitude(C_word k, C_word result)
2353{
2354  if (C_bignum_negativep(result)) {
2355    C_word *scan, *end, digit, carry;
2356
2357    scan = C_bignum_digits(result);
2358    end = scan + C_bignum_size(result);
2359    carry = 1;
2360
2361    while (scan < end) {
2362      digit = (~*scan & BIGNUM_DIGIT_MASK) + carry;
2363
2364      if (C_fitsinbignumdigitp(digit)) {
2365        carry = 0;
2366      } else {
2367        digit &= C_BIGNUM_DIGIT_MASK;
2368        carry = 1;
2369      }
2370   
2371      *scan++ = digit;
2372    }
2373  }
2374  C_kontinue(k, C_bignum_normalize(result));
2375}
2376
2377void C_ccall
2378C_u_bignum_remainder_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y)
2379{
2380  C_word negp = C_mk_bool(C_bignum_negativep(x));
2381   
2382  y = C_unfix(y);
2383
2384  if (y == 1 || y == -1) {
2385    C_kontinue(k, C_fix(0));
2386  } else if (y == C_MOST_NEGATIVE_FIXNUM) {
2387    /* This is the only case we need to go allocate a bignum for */
2388    C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2, size;
2389
2390    y = C_a_u_i_fix_to_big(&a, C_fix(C_MOST_NEGATIVE_FIXNUM));
2391
2392    /* We can skip bignum_divide_2_unsigned because we need no quotient */
2393    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned_2, k, x, y,
2394                   /* Do not return quotient, do return remainder */
2395                   C_SCHEME_FALSE, C_SCHEME_TRUE, negp,
2396                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
2397    size = C_fix(C_bignum_size(x) + 1);
2398    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2399  } else {
2400    C_word absy = (y < 0) ? -y : y;
2401     
2402    if (C_fitsinbignumhalfdigitp(absy)) {
2403      bignum_destructive_remainder_unsigned_halfdigit(k, x, absy, negp);
2404    } else {
2405      C_word kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2,
2406             size = C_fix(C_bignum_size(x)) + 1; /* Due to normalization */
2407     
2408      k2 = C_closure(&ka, 7, (C_word)bignum_destructive_divide_unsigned_digit,
2409                     k, x, C_fix(absy),
2410                     /* Do not return quotient, do return remainder */
2411                     C_SCHEME_FALSE, C_SCHEME_TRUE, negp);
2412      C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2413    }
2414  }
2415}
2416
2417void C_ccall
2418C_u_bignum_remainder_bignum(C_word c, C_word self, C_word k, C_word x, C_word y)
2419{
2420  C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2, size, negp;
2421
2422  switch(bignum_cmp_unsigned(x, y)) {
2423  case 0:
2424    C_kontinue(k, C_fix(0));
2425  case -1:
2426    C_kontinue(k, x);
2427  case 1:
2428    negp = C_bignum_negativep(x);
2429
2430    /* We can skip bignum_divide_2_unsigned because we need no quotient */
2431    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned_2, k, x, y,
2432                   /* Do not return quotient, do return remainder */
2433                   C_SCHEME_FALSE, C_SCHEME_TRUE, negp,
2434                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
2435    size = C_fix(C_bignum_size(x) + 1); /* May need to be normalized */
2436    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2437  }
2438}
2439
2440static void
2441bignum_destructive_remainder_unsigned_halfdigit(C_word k, C_word n, C_word d, C_word negp)
2442{
2443  C_word two_digits,
2444         *start = C_bignum_digits(n),
2445         *scan = start + C_bignum_size(n),
2446         r = 0;
2447
2448  assert((d > 1) && (d < BIGNUM_RADIX_ROOT));
2449  while (start < scan) {
2450    two_digits = (*--scan);
2451    r = C_BIGNUM_DIGIT_COMBINE(r, C_BIGNUM_DIGIT_HI_HALF(two_digits)) % d;
2452    r = C_BIGNUM_DIGIT_COMBINE(r, C_BIGNUM_DIGIT_LO_HALF(two_digits)) % d;
2453  }
2454  C_kontinue(k, C_truep(negp) ? C_fix(-r) : C_fix(r));
2455}
2456
2457void C_ccall
2458C_u_bignum_quotient_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y)
2459{
2460  y = C_unfix(y);
2461 
2462  if (y == 1) {
2463    C_kontinue(k, x);
2464  } else if (y == -1) {
2465    C_u_bignum_negate(1, (C_word)NULL, k, x);
2466  } else if (y == C_MOST_NEGATIVE_FIXNUM) {
2467    /* This is the only case we need to go allocate a bignum for */
2468    C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2,
2469           negp = C_mk_bool(!C_bignum_negativep(x)), size;
2470
2471    y = C_a_u_i_fix_to_big(&a, C_fix(C_MOST_NEGATIVE_FIXNUM));
2472
2473    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned, k, x, y,
2474                   /* Return quotient, not remainder */
2475                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE,
2476                   /* Will be filled in later */
2477                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
2478    size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
2479    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2480  } else {
2481    C_word negp = C_mk_bool((y < 0) ? !C_bignum_negativep(x) :
2482                                      C_bignum_negativep(x)),
2483           absy = (y < 0) ? -y : y,
2484           kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2,
2485           size, func;
2486
2487    if (C_fitsinbignumhalfdigitp(absy)) {
2488      size = C_fix(C_bignum_size(x));
2489      func = (C_word)bignum_destructive_divide_unsigned_halfdigit;
2490    } else {
2491      size = C_fix(C_bignum_size(x)) + 1; /* Due to normalization */
2492      func = (C_word)bignum_destructive_divide_unsigned_digit;
2493    }
2494
2495    k2 = C_closure(&ka, 7, func, k, x, C_fix(absy),
2496                   /* Return quotient, not remainder */
2497                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
2498    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2499  }
2500}
2501
2502void C_ccall
2503C_u_bignum_quotient_bignum(C_word c, C_word self, C_word k, C_word x, C_word y)
2504{
2505  C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2, size,
2506         negp = C_mk_bool(C_bignum_negativep(x) ?
2507                         !C_bignum_negativep(y) :
2508                         C_bignum_negativep(y));
2509
2510  switch(bignum_cmp_unsigned(x, y)) {
2511  case 0:
2512    C_kontinue(k, C_truep(negp) ? C_fix(-1) : C_fix(1));
2513  case -1:
2514    C_kontinue(k, C_fix(0));
2515  case 1:
2516    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned, k, x, y,
2517                   /* Return quotient, not remainder */
2518                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE,
2519                   /* Will be filled in later */
2520                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
2521    size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
2522    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2523  }
2524}
2525
2526static void
2527bignum_destructive_divide_unsigned_halfdigit(C_word c, C_word self, C_word quotient)
2528{
2529  C_word k = C_block_item(self, 1),
2530         numerator = C_block_item(self, 2),
2531         denominator = C_unfix(C_block_item(self, 3)),
2532         /* return_quotient = C_block_item(self, 4), */
2533         return_remainder = C_block_item(self, 5),
2534         remainder_negp = C_block_item(self, 6),
2535         *start = C_bignum_digits(quotient),
2536         *end = start + C_bignum_size(quotient) - 1,
2537         remainder;
2538
2539  bignum_digits_destructive_copy(quotient, numerator);
2540
2541  remainder = bignum_digits_destructive_scale_down(start, end, denominator);
2542  assert(C_fitsinbignumhalfdigitp(remainder));
2543
2544  quotient = C_bignum_normalize(quotient);
2545
2546  if (C_truep(return_remainder)) {
2547    remainder = C_truep(remainder_negp) ? -remainder : remainder;
2548    C_values(4, C_SCHEME_UNDEFINED, k, quotient, C_fix(remainder));
2549  } else {
2550    C_kontinue(k, quotient);
2551  }
2552}
2553
2554/* XXX TODO: This seems unnecessary */
2555static void
2556bignum_destructive_divide_unsigned_digit(C_word c, C_word self, C_word quotient)
2557{
2558  C_word k = C_block_item(self, 1),
2559         numerator = C_block_item(self, 2),
2560         denominator = C_unfix(C_block_item(self, 3)),
2561         return_quotient = C_block_item(self, 4),
2562         return_remainder = C_block_item(self, 5),
2563         remainder_negp = C_block_item(self, 6),
2564         *start = C_bignum_digits(quotient),
2565         *scan = start + C_bignum_size(quotient),
2566         qj, remainder = 0, shift = 0;
2567
2568  /* Because `bignum_divide_digit' requires a normalized denominator. */
2569  while (denominator < ((C_word)1 << (C_BIGNUM_DIGIT_LENGTH-1))) {
2570    denominator <<= 1;
2571    shift++;
2572  }
2573
2574  if (shift) {
2575    bignum_destructive_normalize(quotient, numerator, shift);
2576  } else {
2577    /* Could alo use destructive_normalize here, but this is faster */
2578    bignum_digits_destructive_copy(quotient, numerator);
2579    *(scan-1) = 0; /* Most significant digit, not copied */
2580  }
2581
2582  if (C_truep(return_quotient)) {
2583    while (start < scan) {
2584      remainder = bignum_divide_digit(remainder, (*--scan), denominator, &qj);
2585      (*scan) = qj;
2586    }
2587    quotient = C_bignum_normalize(quotient);
2588  } else {
2589    while (start < scan)
2590      remainder = bignum_divide_digit(remainder, (*--scan), denominator, &qj);
2591  }
2592
2593  if (C_truep(return_remainder)) {
2594    remainder >>= shift;
2595    remainder = C_truep(remainder_negp) ? -remainder : remainder;
2596
2597    if (C_truep(return_quotient))
2598      C_values(4, C_SCHEME_UNDEFINED, k, quotient, C_fix(remainder));
2599    else
2600      C_kontinue(k, C_fix(remainder));
2601  } else {
2602    assert(C_truep(return_quotient));
2603    C_kontinue(k, quotient);
2604  }
2605}
2606
2607static void
2608bignum_destructive_normalize(C_word target, C_word source, C_word shift_left)
2609{
2610  C_word digit, carry = 0,
2611         *scan_source = C_bignum_digits(source),
2612         *scan_target = C_bignum_digits(target),
2613         *end_source = scan_source + C_bignum_size(source),
2614         *end_target = scan_target + C_bignum_size(target),
2615         shift_right = C_BIGNUM_DIGIT_LENGTH - shift_left,
2616         mask = (1L << shift_right) - 1;
2617
2618  while (scan_source < end_source) {
2619    digit = (*scan_source++);
2620    (*scan_target++) = (((digit & mask) << shift_left) | carry);
2621    carry = (digit >> shift_right);
2622  }
2623  if (scan_target < end_target)
2624    (*scan_target) = carry;
2625  else
2626    assert(carry == 0);
2627}
2628
2629/* TODO: This pointer stuff here is suspicious: it's probably slow */
2630static C_word
2631bignum_divide_and_subtract_digit(C_word v1, C_word v2, C_word guess, C_word *u)
2632{
2633  C_word product, diff, sum, carry;
2634
2635  product = v2 * guess;
2636  diff = u[2] - C_BIGNUM_DIGIT_LO_HALF(product);
2637  if (diff < 0) {
2638    u[2] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
2639    carry = C_BIGNUM_DIGIT_HI_HALF(product) + 1;
2640  } else {
2641    u[2] = diff;
2642    carry = C_BIGNUM_DIGIT_HI_HALF(product);
2643  }
2644 
2645  product = v1 * guess + carry;
2646  diff = u[1] - C_BIGNUM_DIGIT_LO_HALF(product);
2647  if (diff < 0) {
2648    u[1] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
2649    carry = C_BIGNUM_DIGIT_HI_HALF(product) + 1;
2650  } else {
2651    u[1] = diff;
2652    carry = C_BIGNUM_DIGIT_HI_HALF(product);
2653    if (carry == 0) return guess; /* DONE */
2654  }
2655
2656  diff = u[0] - carry;
2657  if (diff < 0) {
2658    u[0] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
2659  } else {
2660    u[0] = diff;
2661    return guess; /* DONE */
2662  }
2663
2664  sum = v2 + u[2];
2665  u[2] = sum & C_BIGNUM_HALF_DIGIT_MASK;
2666  carry = C_fitsinbignumhalfdigitp(sum);
2667
2668  sum = v1 + u[1] + carry;
2669  u[1] = sum & C_BIGNUM_HALF_DIGIT_MASK;
2670  carry = C_fitsinbignumhalfdigitp(sum);
2671
2672  u[0] += carry;
2673
2674  return guess - 1;
2675}
2676
2677/* This is a reduced version of the division algorithm, applied to the
2678   case of dividing two bignum digits by one bignum digit.  It is
2679   assumed that the numerator, denominator are normalized. */
2680
2681#define BDD_STEP(qn, j)                                                 \
2682{                                                                       \
2683  uj = u[j];                                                            \
2684  if (uj != v1) {                                                       \
2685    uj_uj1 = C_BIGNUM_DIGIT_COMBINE(uj, u[j + 1]);                      \
2686    guess = uj_uj1 / v1;                                                \
2687    comparand = C_BIGNUM_DIGIT_COMBINE(uj_uj1 % v1, u[j + 2]);          \
2688  } else {                                                              \
2689    guess = C_BIGNUM_HALF_DIGIT_MASK;                                   \
2690    comparand = C_BIGNUM_DIGIT_COMBINE(u[j + 1] + v1, u[j + 2]);        \
2691  }                                                                     \
2692  while ((guess * v2) > comparand) {                                    \
2693    guess -= 1;                                                         \
2694    comparand += v1 << C_BIGNUM_HALF_DIGIT_LENGTH;                      \
2695    if (!C_fitsinbignumdigitp(comparand))                               \
2696      break;                                                            \
2697  }                                                                     \
2698  qn = bignum_divide_and_subtract_digit(v1, v2, guess, &u[j]);          \
2699}
2700
2701/* Because the algorithm from Knuth requires combining the two highest
2702 * digits of u (which we can't fit in a C_word), we instead do two
2703 * such steps, a halfdigit at a time.
2704 */
2705static C_word
2706bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q)
2707{
2708  C_word guess, comparand, v1, v2, uj, uj_uj1, q1, q2, u[4];
2709
2710  if (uh == 0) {
2711    if (ul < v) {
2712      *q = 0;
2713      return ul;
2714    } else if (ul == v) {
2715      *q = 1;
2716      return 0;
2717    }
2718  }
2719  u[0] = C_BIGNUM_DIGIT_HI_HALF(uh);
2720  u[1] = C_BIGNUM_DIGIT_LO_HALF(uh);
2721  u[2] = C_BIGNUM_DIGIT_HI_HALF(ul);
2722  u[3] = C_BIGNUM_DIGIT_LO_HALF(ul);
2723  v1 = C_BIGNUM_DIGIT_HI_HALF(v);
2724  v2 = C_BIGNUM_DIGIT_LO_HALF(v);
2725  BDD_STEP(q1, 0);
2726  BDD_STEP(q2, 1);
2727  *q = C_BIGNUM_DIGIT_COMBINE(q1, q2);
2728  return C_BIGNUM_DIGIT_COMBINE(u[2], u[3]);
2729}
2730
2731#undef BDD_STEP
2732
2733/* Full bignum division */
2734
2735static void
2736bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient)
2737{
2738  C_word numerator = C_block_item(self, 2),
2739         remainder_negp = C_block_item(self, 6),
2740         size = C_fix(C_bignum_size(numerator) + 1);
2741
2742  /* Nice: We can recycle the current closure */
2743  C_set_block_item(self, 0, (C_word)bignum_divide_2_unsigned_2);
2744  C_set_block_item(self, 7, quotient);
2745  C_allocate_bignum(3, (C_word)NULL, self, size, remainder_negp, C_SCHEME_FALSE);
2746}
2747
2748/* For help understanding this algorithm, see:
2749   Knuth, Donald E., "The Art of Computer Programming",
2750   volume 2, "Seminumerical Algorithms"
2751   section 4.3.1, "Multiple-Precision Arithmetic". */
2752
2753static void
2754bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder)
2755{
2756  C_word k = C_block_item(self, 1),
2757         numerator = C_block_item(self, 2),
2758         denominator = C_block_item(self, 3),
2759         return_quotient = C_block_item(self, 4),
2760         return_remainder = C_block_item(self, 5),
2761         /* This one may be overwritten with the remainder */
2762         /* remainder_negp = C_block_item(self, 6), */
2763         quotient = C_block_item(self, 7),
2764         length_d = C_bignum_size(denominator),
2765         d1 = *(C_bignum_digits(denominator) + length_d - 1),
2766         shift = 0;
2767
2768  assert(length_d > 1);
2769  while (d1 < (BIGNUM_RADIX / 2)) {
2770    d1 <<= 1;
2771    shift++;
2772  }
2773
2774  if (shift == 0) { /* Already normalized */
2775    bignum_digits_destructive_copy(remainder, numerator);
2776    /* Set most significant digit */
2777    *(C_bignum_digits(remainder) + C_bignum_size(numerator)) = 0;
2778
2779    bignum_destructive_divide_normalized(remainder, denominator, quotient);
2780
2781    if (C_truep(C_and(return_quotient, return_remainder))) {
2782      C_values(4, C_SCHEME_UNDEFINED, k,
2783               C_bignum_normalize(quotient), C_bignum_normalize(remainder));
2784    } else if (C_truep(return_remainder)) {
2785      C_kontinue(k, C_bignum_normalize(remainder));
2786    } else {
2787      assert(C_truep(return_quotient));
2788      C_kontinue(k, C_bignum_normalize(quotient));
2789    }
2790  } else {
2791    /* Requires normalisation of denominator;  Allocate temp bignum for it. */
2792    C_set_block_item(self, 0, (C_word)bignum_divide_2_unsigned_3);
2793    C_set_block_item(self, 6, remainder);
2794    C_set_block_item(self, 8, C_fix(shift));
2795    C_allocate_bignum(3, (C_word)NULL, self, C_fix(length_d),
2796                      C_SCHEME_FALSE, C_SCHEME_FALSE);
2797  }
2798}
2799
2800static void
2801bignum_divide_2_unsigned_3(C_word c, C_word self, C_word tmp_big)
2802{
2803  C_word k = C_block_item(self, 1),
2804         numerator = C_block_item(self, 2),
2805         denominator = C_block_item(self, 3),
2806         return_quotient = C_block_item(self, 4),
2807         return_remainder = C_block_item(self, 5),
2808         remainder = C_block_item(self, 6),
2809         quotient = C_block_item(self, 7),
2810         shift = C_unfix(C_block_item(self, 8));
2811
2812  bignum_destructive_normalize(remainder, numerator, shift);
2813  bignum_destructive_normalize(tmp_big, denominator, shift);
2814  bignum_destructive_divide_normalized(remainder, tmp_big, quotient);
2815
2816  if (C_truep(return_remainder)) {
2817    if (C_truep(return_quotient)) {
2818      C_values(4, C_SCHEME_UNDEFINED, k,
2819               C_bignum_normalize(quotient),
2820               bignum_normalize_shifted(remainder, shift));
2821    } else {
2822      C_kontinue(k, bignum_normalize_shifted(remainder, shift));
2823    }
2824  } else {
2825    assert(C_truep(return_quotient));
2826    C_kontinue(k, C_bignum_normalize(quotient));
2827  }
2828}
2829
2830static void
2831bignum_destructive_divide_normalized(C_word u, C_word v, C_word q)
2832{
2833  C_word u_length = C_bignum_size(u),
2834         v_length = C_bignum_size(v),
2835         *u_start = C_bignum_digits(u),
2836         *u_scan = u_start + u_length,
2837         *u_scan_limit = u_start + v_length,
2838         *u_scan_start = u_scan - v_length,
2839         *v_start = C_bignum_digits(v),
2840         *v_end = v_start + v_length,
2841         *q_scan = (q == C_SCHEME_UNDEFINED) ? NULL :
2842                   (C_bignum_digits(q) + C_bignum_size(q)),
2843         v1 = v_end[-1],
2844         v2 = v_end[-2],
2845         ph, /* high half of double-digit product */
2846         pl, /* low half of double-digit product */
2847         guess, uj, qj,
2848         gh, /* high half-digit of guess */
2849         ch, /* high half of double-digit comparand */
2850         v2l = C_BIGNUM_DIGIT_LO_HALF(v2),
2851         v2h = C_BIGNUM_DIGIT_HI_HALF(v2),
2852         cl, /* low half of double-digit comparand */
2853         gl, /* low half-digit of guess */
2854         gm; /* memory loc for reference parameter */
2855
2856  while (u_scan_limit < u_scan) {
2857    uj = (*--u_scan);
2858    if (uj != v1) {
2859      /* comparand = (((((uj * B) + uj1) % v1) * b) + uj2);
2860         guess = (((uj * B) + uj1) / v1); */
2861      cl = u_scan[-2];
2862      ch = bignum_divide_digit(uj, (u_scan[-1]), v1, (&gm));
2863      guess = gm;
2864    } else {
2865      cl = u_scan[-2];
2866      ch = u_scan[-1] + v1;
2867      guess = C_BIGNUM_DIGIT_MASK;
2868    }
2869    while (1) {
2870      /* product = (guess * v2); */
2871      gl = C_BIGNUM_DIGIT_LO_HALF(guess);
2872      gh = C_BIGNUM_DIGIT_HI_HALF(guess);
2873      pl = v2l * gl;
2874      ph = v2l * gh + v2h * gl + C_BIGNUM_DIGIT_HI_HALF(pl);
2875      pl = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(ph),
2876                                  C_BIGNUM_DIGIT_LO_HALF(pl));
2877      ph = v2h * gh + C_BIGNUM_DIGIT_HI_HALF(ph);
2878      /* if (comparand >= product) */
2879      if ((ch > ph) || ((ch == ph) && (cl >= pl)))
2880        break;
2881      guess--;
2882      /* comparand += (v1 << C_BIGNUM_DIGIT_LENGTH) */
2883      ch += v1;
2884      /* if (comparand >= (B^2)) */
2885      if (!C_fitsinbignumdigitp(ch))
2886        break;
2887    }
2888    qj = bignum_divide_and_subtract(v_start, v_end, guess, (--u_scan_start));
2889    if (q_scan != NULL)
2890      (*--q_scan) = qj;
2891  }
2892}
2893
2894static C_word
2895bignum_divide_and_subtract(C_word *v_start, C_word *v_end, C_word guess, C_word *u_start)
2896{
2897  if (guess == 0) {
2898    return 0;
2899  } else {
2900    C_word *v_scan = v_start,
2901           *u_scan = u_start,
2902           carry = 0,
2903           gl, gh, v, pl, vl, vh, ph, diff, sum;
2904
2905    gl = C_BIGNUM_DIGIT_LO_HALF(guess);
2906    gh = C_BIGNUM_DIGIT_HI_HALF(guess);
2907    while (v_scan < v_end) {
2908      v = (*v_scan++);
2909      vl = C_BIGNUM_DIGIT_LO_HALF(v);
2910      vh = C_BIGNUM_DIGIT_HI_HALF(v);
2911      pl = vl * gl + C_BIGNUM_DIGIT_LO_HALF(carry);
2912      ph = vl * gh + vh * gl + C_BIGNUM_DIGIT_HI_HALF(pl) +
2913                               C_BIGNUM_DIGIT_HI_HALF(carry);
2914      diff = (*u_scan) - C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(ph),
2915                                                C_BIGNUM_DIGIT_LO_HALF(pl));
2916      if (diff < 0) {
2917        (*u_scan++) = diff + ((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
2918        carry = vh * gh + C_BIGNUM_DIGIT_HI_HALF(ph) + 1;
2919      } else {
2920        (*u_scan++) = diff;
2921        carry = vh * gh + C_BIGNUM_DIGIT_HI_HALF(ph);
2922      }
2923    }
2924    if (carry == 0)
2925      return guess;
2926
2927    diff = ((*u_scan) - carry);
2928    if (diff < 0) {
2929      (*u_scan) = diff + ((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
2930    } else {
2931      (*u_scan) = diff;
2932      return guess;
2933    }
2934 
2935    /* Subtraction generated carry, implying guess is one too large.
2936       Add v back in to bring it back down. */
2937    v_scan = v_start;
2938    u_scan = u_start;
2939    carry = 0;
2940    while (v_scan < v_end) {
2941      sum = ((*v_scan++) + (*u_scan) + carry);
2942      if (sum < BIGNUM_RADIX) {
2943        (*u_scan++) = sum;
2944        carry = 0;
2945      } else {
2946        (*u_scan++) = sum & C_BIGNUM_DIGIT_MASK;
2947        carry = 1;
2948      }
2949    }
2950    if (carry) {
2951      sum = (*u_scan) + carry;
2952      (*u_scan) = C_fitsinbignumdigitp(sum) ? sum : (sum  & C_BIGNUM_DIGIT_MASK);
2953    }
2954    return guess - 1;
2955  }
2956}
2957
2958/* Like bignum_normalize, but this also shifts division-normalized numbers */
2959static C_word
2960bignum_normalize_shifted(C_word big, C_word shift_right)
2961{
2962  C_word length = C_bignum_size(big),
2963        *start = C_bignum_digits(big),
2964        *scan = start + length,
2965        digit, carry = 0,
2966        shift_left = C_BIGNUM_DIGIT_LENGTH - shift_right,
2967        mask = (1L << shift_right) - 1;
2968 
2969  while (!(*--scan)) {
2970    if (start == scan) { /* Don't bother with anything else */
2971      return C_fix(0);
2972    }
2973    --length;
2974  }
2975
2976  digit = (*scan);
2977  (*scan) = (digit >> shift_right);
2978  length -= (*scan == 0); /* Add 1 or 0 */
2979  carry = ((digit & mask) << shift_left);
2980 
2981  while (start < scan) {
2982    digit = (*--scan);
2983    (*scan) = ((digit >> shift_right) | carry);
2984    carry = ((digit & mask) << shift_left);
2985  }
2986  assert(carry == 0);
2987  assert(C_bignum_size(big) != length);
2988  assert(length != 1 || *C_bignum_digits(big) != 0);
2989
2990  switch(length) {
2991  case 0:
2992    return C_fix(0);
2993  case 1:
2994    return C_fix(C_bignum_negativep(big) ? -*start : *start);
2995  case 2:
2996    if (C_bignum_negativep(big) && *scan == 1 && *start == 0)
2997      return C_fix(C_MOST_NEGATIVE_FIXNUM);
2998    /* FALLTHROUGH */
2999  default:
3000    /* Mutate vector size of internal bignum vector. */
3001    C_block_header(C_internal_bignum(big)) = (C_STRING_TYPE | C_wordstobytes(length+1));
3002    /* Set internal header. */
3003    C_bignum_header(big) = (C_bignum_header(big) & C_BIGNUM_HEADER_SIGN_BIT) | length;
3004    return big;
3005  }
3006}
Note: See TracBrowser for help on using the repository browser.