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

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

numbers: Replace the silly memcpy calls with a function to abstract out the mess

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