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

Last change on this file since 31407 was 31407, checked in by sjamaan, 5 years ago

numbers: Roll trimming and normalizing of bignums into one function, simplifying the returning of bignums, as well as offering a tiny performance improvement.

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