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

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

numbers: Minor simplification by cutting out the ability to turn off assertion checks specifically for bignums. Remove ancient constants BIG_FREE and FORCE_FINALIZERS (which go all the way back to the gmp implementation)

File size: 71.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
42static void *tags;
43
44#define fix_to_flo(p, n, f)       C_flonum(p, C_unfix(f))
45#define flonump(x)                C_mk_bool(!C_immediatep(x) && C_block_header(x) == C_FLONUM_TAG)
46#define big_of(v)                 ((bignum_type)C_data_pointer(C_block_item(v, 1)))
47#define big_negp(b)               C_mk_bool(BIGNUM_NEGATIVE_P(big_of(b)))
48#define big_oddp(b)               C_mk_bool(BIGNUM_REF(big_of(b), 0) & 1)
49
50static C_word
51init_tags(___scheme_value tagvec)
52{
53  tags = CHICKEN_new_gc_root();
54  CHICKEN_gc_root_set(tags, tagvec);
55  return C_SCHEME_UNDEFINED;
56}
57
58/*
59 * This is an odd one out. It doesn't accept a continuation.
60 * I've put in the (more or less) verbatim code for s48_bignum_to_double.
61 */
62static C_word
63big_to_flo(C_word **p, C_word n, C_word value)
64{
65  bignum_type bignum = big_of(value);
66  double accumulator = 0;
67  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
68  bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
69  while (start < scan)
70    accumulator = ((accumulator * BIGNUM_RADIX) + (*--scan));
71  return C_flonum(p, (BIGNUM_NEGATIVE_P(bignum) ? -accumulator : accumulator));
72}
73
74
75static C_word
76check_number(C_word x)
77{
78  C_word tagvec;
79
80  if((x & C_FIXNUM_BIT) != 0) return C_fix(FIX);
81  else if(C_immediatep(x)) return C_fix(NONE);
82  else if(C_header_bits(x) == C_FLONUM_TYPE) return C_fix(FLO);
83  else if(C_header_bits(x) == C_STRUCTURE_TYPE) {
84    tagvec = CHICKEN_gc_root_ref(tags);
85   
86    if (C_block_item(x, 0) == C_block_item(tagvec, BIG_TAG))
87      return C_fix(BIG);
88    else if (C_block_item(x, 0) == C_block_item(tagvec, COMP_TAG))
89      return C_fix(COMP);
90    else if (C_block_item(x, 0) == C_block_item(tagvec, RAT_TAG))
91      return C_fix(RAT);
92    else
93      return C_fix(NONE);
94  } else
95      return C_fix(NONE);
96}
97
98static void
99C_wrap_bignum(C_word c, C_word closure, C_word bigvec)
100{
101  C_word ab[3], *a = ab;
102  bignum_type src = (bignum_type)C_block_item(C_block_item(closure, 2), 0);
103  bignum_digit_type *scan = BIGNUM_TO_POINTER(src);
104  bignum_digit_type *end = BIGNUM_START_PTR(src) + BIGNUM_LENGTH(src);
105  bignum_digit_type *tgt = C_data_pointer(bigvec);
106  C_word tagvec = CHICKEN_gc_root_ref(tags);
107  C_word result = C_structure(&a, 2, C_block_item(tagvec, BIG_TAG), bigvec);
108  while (scan < end)
109    (*tgt++) = (*scan++);
110  BIGNUM_DEALLOCATE(src);
111  C_kontinue(C_block_item(closure, 1), result);
112}
113
114static void
115C_return_bignum(C_word k, bignum_type b)
116{
117  C_word result;
118 
119  assert(b != BIGNUM_OUT_OF_BAND);
120
121  switch(BIGNUM_LENGTH(b)) {
122  case 0:
123    BIGNUM_DEALLOCATE(b);
124    C_kontinue(k, C_fix(0));
125  /* This code "knows" that bignums have 2 "reserved" bits, like fixnums */
126  case 1:
127    result = C_fix(BIGNUM_NEGATIVE_P(b) ? -BIGNUM_REF(b, 0) : BIGNUM_REF(b, 0));
128    BIGNUM_DEALLOCATE(b);
129    C_kontinue(k, result);
130  case 2:
131    /* Edge case: most negative fixnum */
132    if (BIGNUM_NEGATIVE_P(b) && BIGNUM_REF(b, 0) == 0 && BIGNUM_REF(b, 1) == 1) {
133      BIGNUM_DEALLOCATE(b);
134      C_kontinue(k, C_fix(C_MOST_NEGATIVE_FIXNUM));
135    }
136    /* FALLTHROUGH */
137  default:
138    {
139      C_word pab[2], *pa = pab, kab[4], *ka = kab, k2;
140      /* Make result a wrapped pointer because C_closure wants scheme objects */
141      result = C_mpointer(&pa, b);
142      k2 = C_closure(&ka, 3, (C_word)C_wrap_bignum, k, result);
143      /* Here we assume bignum digits are C words.. */
144      C_allocate_vector(6, (C_word)NULL, k2,
145                        C_fix(sizeof(C_word) * (BIGNUM_LENGTH(b) + 1)),
146                        /* Byte vec, no initialization, align at 8 bytes */
147                        C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
148    }
149  }
150}
151
152static void
153C_bignum_wrapped_return_bigobj(C_word c, C_word closure, C_word wrapped_big)
154{
155  C_word obj = C_block_item(closure, 2);
156  C_values(4, C_SCHEME_UNDEFINED, C_block_item(closure, 1), wrapped_big, obj);
157}
158
159static void
160C_return_big_fix(C_word k, bignum_type big, C_word fix)
161{
162  C_word kab[4], *ka = kab, k2;
163  k2 = C_closure(&ka, 3, (C_word)C_bignum_wrapped_return_bigobj, k, fix);
164  C_return_bignum(k2, big);
165}
166
167static void
168C_b2_wrapped(C_word c, C_word closure, C_word wrapped_b2)
169{
170  C_word k = C_block_item(closure, 1); /* Original closure */
171  C_word kab[4], *ka = kab, k2; /* Next continuation */
172  k2 = C_closure(&ka, 3, (C_word)C_bignum_wrapped_return_bigobj, k, wrapped_b2);
173  C_return_bignum(k2, (bignum_type)C_block_item(C_block_item(closure, 2), 0));
174}
175
176static void
177C_return_2_bignums(C_word k, bignum_type b1, bignum_type b2)
178{
179  C_word bab[2], *ba = bab, kab[4], *ka = kab, k2, b1_ptr;
180  /* Make b1 a wrapped pointer because C_closure wants scheme objects */
181  b1_ptr = C_mpointer(&ba, b1);
182  /* Wrap b2 first, then b1. Return them to k in the original (b1,b2) order */
183  k2 = C_closure(&ka, 3, (C_word)C_b2_wrapped, k, b1_ptr);
184  C_return_bignum(k2, b2);
185}
186
187
188static bignum_type
189bignum_allocate(bignum_length_type length, int negative_p)
190{
191  bignum_type result;
192  bignum_digit_type *digits;
193
194  digits = (bignum_digit_type *)C_malloc(sizeof(bignum_digit_type)*(length + 1));
195 
196  if(digits == NULL) {
197    fprintf(stderr, "out of memory - can not allocate number");
198    exit(EXIT_FAILURE);
199  }
200
201  result = (bignum_type)digits;
202  BIGNUM_SET_HEADER(result, length, negative_p);
203  return result;
204}
205
206static bignum_type
207bignum_allocate_zeroed(bignum_length_type length, int negative_p)
208{
209  assert((length >= 0) || (length < BIGNUM_RADIX));
210  {
211    bignum_type result;
212    bignum_digit_type *digits;
213    bignum_digit_type *scan;
214    bignum_digit_type *end;
215    digits=(bignum_digit_type *)C_malloc(sizeof(bignum_digit_type)*(length+1));
216 
217    if(digits == NULL) {
218      fprintf(stderr, "out of memory - can not allocate number");
219      exit(EXIT_FAILURE);
220    }
221
222    result = (bignum_type)digits;
223    BIGNUM_SET_HEADER(result, length, negative_p);
224    scan = BIGNUM_START_PTR(result);
225    end = scan + length;
226    while (scan < end)
227      (*scan++) = 0;
228    return (result);
229  }
230}
231
232static bignum_type
233bignum_allocate_from_fixnum(C_word fix)
234{
235  bignum_type ret;
236
237  if (fix == C_fix(0)) {
238    ret = bignum_allocate(0, 0);
239  } else if (fix == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
240    ret = bignum_allocate(2, 1);
241    BIGNUM_REF(ret, 0) = 0;
242    BIGNUM_REF(ret, 1) = 1;
243  } else {
244    bignum_digit_type digit = C_unfix(fix);
245    ret = bignum_allocate(1, digit < 0);
246    BIGNUM_REF(ret, 0) = ((digit < 0) ? -digit : digit);
247  }
248  return ret;
249}
250
251static bignum_type
252bignum_digit_to_bignum(bignum_digit_type digit, int neg_p)
253{
254  bignum_type ret;
255  if (digit == 0)
256    return bignum_allocate(0, 0);
257
258  ret = bignum_allocate(1, neg_p);
259  BIGNUM_REF(ret, 0) = digit;
260  return ret;
261}
262
263static bignum_type
264shorten_bignum(bignum_type big, bignum_length_type newlength)
265{
266  bignum_digit_type *digits, *newdigits;
267  digits = BIGNUM_TO_POINTER(big);
268  newdigits = (bignum_digit_type *)C_realloc(digits, sizeof(bignum_digit_type)*(newlength + 1));
269  if (newdigits == NULL) {
270    fprintf(stderr, "out of memory - can not reallocate number");
271    exit(EXIT_FAILURE);
272  }
273  return (bignum_type)newdigits;
274}
275
276static bignum_type
277bignum_trim(bignum_type bignum)
278{
279  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
280  bignum_digit_type * end = (start + (BIGNUM_LENGTH (bignum)));
281  bignum_digit_type * scan = end;
282  while ((start <= scan) && ((*--scan) == 0))
283    ;
284  scan += 1;
285  if (scan < end)
286    {
287      bignum_length_type length = (scan - start);
288      BIGNUM_SET_HEADER
289        (bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
290      BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
291    }
292  return (bignum);
293}
294
295/* Copying */
296
297static bignum_type
298bignum_copy(bignum_type source)
299{
300  bignum_type target =
301    (bignum_allocate ((BIGNUM_LENGTH (source)), (BIGNUM_NEGATIVE_P (source))));
302  bignum_destructive_copy (source, target);
303  return (target);
304}
305
306static void
307bignum_destructive_copy(bignum_type source, bignum_type target)
308{
309  bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
310  bignum_digit_type * end_source =
311    (scan_source + (BIGNUM_LENGTH (source)));
312  bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
313  while (scan_source < end_source)
314    (*scan_target++) = (*scan_source++);
315  return;
316}
317
318static bignum_type
319bignum_new_sign(bignum_type bignum, int negative_p)
320{
321  bignum_type result =
322    (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
323  bignum_destructive_copy (bignum, result);
324  return (result);
325}
326
327#define fix_randomize  C_randomize
328#define fix_random C_random_fixnum
329
330/*
331 * This random number generator is very simple. Probably too simple...
332 */
333static C_word
334big_randomize(C_word bignum)
335{
336  bignum_type big;
337  C_word seed = 0;
338  bignum_digit_type * scan, * end;
339
340  big = big_of(bignum);
341  scan = (BIGNUM_START_PTR (big));
342  end = scan + BIGNUM_LENGTH(big);
343  /* What a cheap way to initialize the random generator. I feel dirty! */
344  while (scan < end)
345    seed ^= *scan++;
346
347  srand(C_unfix(seed));
348  return C_SCHEME_UNDEFINED;
349}
350
351static void
352big_random(C_word c, C_word self, C_word k, C_word max)
353{
354  bignum_type bigmax;
355  bignum_type result;
356  bignum_length_type randlen, max_len, max_bits;
357  bignum_digit_type max_top_digit, d;
358  bignum_digit_type * scan, * end;
359 
360  bigmax = big_of(max);
361
362  max_len = BIGNUM_LENGTH(bigmax);
363  max_top_digit = d = BIGNUM_REF(bigmax, max_len - 1);
364 
365  max_bits = (max_len - 1) * BIGNUM_DIGIT_LENGTH;
366  while(d) {
367    max_bits++;
368    d >>= 1;
369  }
370  /* Subtract/add one because we don't want zero to be over-represented */
371  randlen = (((double)rand())/(RAND_MAX + 1.0) * (double)(max_bits - 1));
372  randlen = BIGNUM_BITS_TO_DIGITS(randlen + 1);
373 
374  result = bignum_allocate(randlen, BIGNUM_NEGATIVE_P(bigmax));
375  scan = BIGNUM_START_PTR(result);
376  end = scan + randlen - 1; /* Go to just before the end */
377  while(scan < end)
378    *scan++ = (((double)rand())/(RAND_MAX + 1.0) * (double)BIGNUM_RADIX);
379  /*
380   * Last word is special when length is max_len: It must be less than
381   * max's most significant digit, instead of BIGNUM_RADIX.
382   */
383  if (max_len == randlen)
384    *scan = (((double)rand())/(RAND_MAX + 1.0) * (double)max_top_digit);
385  else
386    *scan = (((double)rand())/(RAND_MAX + 1.0) * (double)BIGNUM_RADIX);
387
388  result = bignum_trim(result); 
389  C_return_bignum(k, result);
390}
391
392static void
393fix_plus_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
394{
395  C_word z;
396
397  /* Exceptional situation: this will cause a real overflow */
398  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) && y == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
399    bignum_type big;
400    bignum_digit_type *ref;
401    big = bignum_allocate(2, 1);
402    ref = BIGNUM_START_PTR(big);
403    *ref++ = 0;
404    *ref = 2;
405    C_return_bignum(k, big);
406  }
407
408  z = C_unfix(x) + C_unfix(y);
409
410  /* This code "knows" that both fixnums and bignums have 2 reserved bits */
411  if(!C_fitsinfixnump(z)) {
412    bignum_type big;
413    bignum_digit_type *ref;
414    big = bignum_allocate(2, (z < 0));
415    ref = BIGNUM_START_PTR(big);
416    *ref++ = labs(z) & BIGNUM_DIGIT_MASK;
417    *ref = 1;
418    C_return_bignum(k, big);
419  }
420
421  C_kontinue(k, C_fix(z));
422}
423
424
425static void
426fix_plus_big(C_word c, C_word self, C_word k, C_word x, C_word y)
427{
428  bignum_type bigx, result;
429  bigx = bignum_allocate_from_fixnum(x);
430  result = bignum_add(bigx, big_of(y));
431  BIGNUM_DEALLOCATE(bigx);
432  C_return_bignum(k, result);
433}
434
435
436static void
437fix_minus_big(C_word c, C_word self, C_word k, C_word x, C_word y)
438{
439  bignum_type tmp, result;
440  tmp = bignum_allocate_from_fixnum(x);
441
442  result = bignum_subtract(tmp, big_of(y));
443  BIGNUM_DEALLOCATE(tmp);
444  C_return_bignum(k, result);
445}
446
447
448static void
449big_minus_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
450{
451  bignum_type tmp, result;
452  tmp = bignum_allocate_from_fixnum(y);
453
454  result = bignum_subtract(big_of(x), tmp);
455  BIGNUM_DEALLOCATE(tmp);
456  C_return_bignum(k, result);
457}
458
459static void
460big_divrem_big(C_word c, C_word self, C_word k, C_word x, C_word y)
461{
462  bignum_type numerator = big_of(x);
463  bignum_type denominator = big_of(y);
464  int q_neg_p = ((BIGNUM_NEGATIVE_P (denominator))
465                 ? (! (BIGNUM_NEGATIVE_P (numerator)))
466                 : (BIGNUM_NEGATIVE_P (numerator)));
467
468  switch (bignum_compare_unsigned (numerator, denominator))
469    {
470    case bignum_comparison_equal:
471      C_values(4, C_SCHEME_UNDEFINED, k,
472               q_neg_p ? C_fix(-1) : C_fix(1), C_fix(0));
473    case bignum_comparison_less:
474      C_values(4, C_SCHEME_UNDEFINED, k, C_fix(0), x);
475    case bignum_comparison_greater:
476    default:                                  /* to appease gcc -Wall */
477      {
478        bignum_type quotient, remainder;
479        int r_neg_p = BIGNUM_NEGATIVE_P (numerator) ? 1 : 0;
480       
481        assert(BIGNUM_LENGTH(denominator) > 1);
482        bignum_divide_unsigned_large_denominator
483          (numerator, denominator, (&quotient), (&remainder), q_neg_p, r_neg_p);
484        C_return_2_bignums(k, quotient, remainder);
485      }
486    }
487}
488
489static void
490big_divrem_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
491{
492  bignum_type bigx = big_of(x);
493  y = C_unfix(y);
494
495  if (y == 1)
496    C_values(4, C_SCHEME_UNDEFINED, k, x, C_fix(0));
497  else if (y == -1)
498    C_return_big_fix(k, bignum_new_sign(bigx, !(BIGNUM_NEGATIVE_P(bigx))), C_fix(0));
499
500  /* Too bad, we really need to do some work... */
501  {
502    int q_neg_p = (y < 0) ? !(BIGNUM_NEGATIVE_P(bigx)) : BIGNUM_NEGATIVE_P(bigx);
503    int r_neg_p = BIGNUM_NEGATIVE_P(bigx);
504    bignum_digit_type abs_y = (y < 0) ? -y : y;
505    bignum_type quotient, remainder;
506   
507    if (y == C_MOST_NEGATIVE_FIXNUM) {
508      if (!BIGNUM_NEGATIVE_P(bigx) && BIGNUM_LENGTH(bigx) == 1
509          && BIGNUM_REF(bigx, 1) == 1 && BIGNUM_REF(bigx, 0) == 0) {
510        /*
511         * Very very special case:
512         * quotient(MOST_NEGATIVE_FIXNUM, -(MOST_NEGATIVE_FIXNUM)) => -1
513         */
514        C_values(4, C_SCHEME_UNDEFINED, k, C_fix(-1), C_fix(0));
515      } else {
516        /* This is the only case we need to go allocate a bignum for */
517        bignum_type bigy =
518          bignum_allocate_from_fixnum(C_fix(C_MOST_NEGATIVE_FIXNUM));
519
520        bignum_divide_unsigned_large_denominator
521          (bigx, bigy, (&quotient), (&remainder), q_neg_p, r_neg_p);
522        BIGNUM_DEALLOCATE(bigy);
523        C_return_2_bignums(k, quotient, remainder);
524      }
525    } else if (abs_y < BIGNUM_RADIX_ROOT) {
526      bignum_divide_unsigned_small_denominator
527        (bigx, abs_y, (&quotient), (&remainder), q_neg_p, r_neg_p);
528      C_return_2_bignums(k, quotient, remainder);
529    } else {
530      bignum_divide_unsigned_medium_denominator
531        (bigx, abs_y, (&quotient), (&remainder), q_neg_p, r_neg_p);
532      C_return_2_bignums(k, quotient, remainder);
533    }
534  }
535}
536
537static void
538C_fixnum_gcd(C_word c, C_word self, C_word k, C_word x, C_word y)
539{
540   C_word r;
541   
542   x = C_unfix(x);
543   y = C_unfix(y);
544   
545   if (x < 0) x = -x;
546   if (y < 0) y = -y;
547   
548   while(y != 0) {
549     r = x % y;
550     x = y;
551     y = r;
552   }
553   C_kontinue(k, C_fix(x));
554}
555
556/*
557 * Big_gcd is a huge function and it sucks that it needs to be in C.
558 *
559 * Why does it need to be in C?  Because if you have a very big bignum
560 * that doesn't cleanly divide another big bignum, you end up calling
561 * the remainder procedure a lot in Scheme.  This produces tons of
562 * intermediate bignums, which means a lot of copies into GC'able memory
563 * need to be made (and the GC will be triggered more often).
564 * That's a major slowdown.  Doing the loop in C means the intermediate
565 * results can be cleaned up right away each loop step, and returning
566 * just one result to Scheme.
567 * Once (if?) we find a way to avoid copying bignums (instead allocating
568 * directly in GCable memory) this function can be cut out and replaced by
569 * a recursive call to gcd-0 in Scheme.
570 */
571static void
572big_gcd(C_word c, C_word self, C_word k, C_word x, C_word y)
573{
574  bignum_type bigx = big_of(x), bigy = big_of(y), bigr;
575
576  assert(BIGNUM_LENGTH(bigx) > 1);
577  assert(BIGNUM_LENGTH(bigy) > 1);
578 
579  switch(bignum_compare_unsigned (bigx, bigy)) {
580  case bignum_comparison_equal:
581    C_kontinue(k, x);
582  case bignum_comparison_less:
583    /* Swap, since remainder of bigx, bigy would be bigx, causing an extra loop */
584    bigr = bigy;
585    bigy = bigx;
586    bigx = bigr;
587
588    /* Make x and y match for the special case where gcd(x, y) = y */
589    {
590      C_word tmp = y;
591      y = x;
592      x = tmp;
593    }
594    /* FALLTHROUGH */
595  default: /* Continue below */
596    break;
597  }
598 
599  /*
600   * Be careful! Don't deallocate live objects. We could start with a copy
601   * or compare pointers with big_of(x) or y every time but that seems wasteful.
602   */
603  bignum_divide_unsigned_large_denominator
604   (bigx, bigy, ((bignum_type *) 0), (&bigr), 0, 0);
605  bigx = bigy;
606  bigy = bigr;
607  /* Original bigx is forgotten now */
608  assert(bigx != big_of(x));
609  assert(bigy != big_of(x));
610  /* Only bigx points to y */
611  assert(bigy != big_of(y));
612  assert(bigx == big_of(y));
613
614  switch (BIGNUM_LENGTH(bigy)) {
615  case 0:  /* bigy = 0 */
616    /* remainder(x, y) = 0 => y  |  x < y */
617    BIGNUM_DEALLOCATE(bigy); /* Got allocated in previous step */
618    C_kontinue(k, y);
619  case 1:
620    if (BIGNUM_REF(bigy, 0) == 1) { /* y = 1? Then 1 is the result */
621      BIGNUM_DEALLOCATE(bigy); /* Got allocated in previous step */
622      C_kontinue(k, C_fix(1));
623    } else if (BIGNUM_REF(bigy, 0) < BIGNUM_RADIX_ROOT)
624      bigr = bignum_remainder_unsigned_small_denominator
625              (bigx, BIGNUM_REF(bigy, 0), 0);
626    else
627      bignum_divide_unsigned_medium_denominator
628       (bigx, BIGNUM_REF(bigy, 0), (bignum_type *)0, (&bigr), 0, 0);
629    break;
630  default:
631    bignum_divide_unsigned_large_denominator
632     (bigx, bigy, ((bignum_type *) 0), (&bigr), 0, 0);
633  }
634  /* Swap, but don't deallocate x since it holds the original value of y */
635  bigx = bigy;
636  bigy = bigr;
637 
638  /* Original bigy is forgotten now, we can safely always deallocate bigx */
639
640  /* Assume that bignums coming from outside are never length 1 */
641  assert(bigx != big_of(y));
642 
643  while(BIGNUM_LENGTH(bigy) > 1) {
644    bignum_divide_unsigned_large_denominator
645     (bigx, bigy, ((bignum_type *) 0), (&bigr), 0, 0);
646    BIGNUM_DEALLOCATE(bigx);
647    bigx = bigy;
648    bigy = bigr;
649  }
650
651  /* Finish up with a faster loop until y = 0 (ie, length(bigy) = 0) */
652  while (BIGNUM_LENGTH(bigy) == 1) {
653    if (BIGNUM_REF(bigy, 0) == 1) {
654      BIGNUM_DEALLOCATE(bigx);
655      BIGNUM_DEALLOCATE(bigy);
656      C_kontinue(k, C_fix(1));
657      break;
658    } else if (BIGNUM_REF(bigy, 0) < BIGNUM_RADIX_ROOT) {
659      bigr = bignum_remainder_unsigned_small_denominator
660              (bigx, BIGNUM_REF(bigy, 0), 0);
661    } else {
662      bignum_divide_unsigned_medium_denominator
663       (bigx, BIGNUM_REF(bigy, 0), (bignum_type *)0, (&bigr), 0, 0);
664    }
665    BIGNUM_DEALLOCATE(bigx);
666    bigx = bigy;
667    bigy = bigr;
668  }
669  BIGNUM_DEALLOCATE(bigy);
670  C_return_bignum(k, bigx);
671}
672
673static C_word
674C_a_i_flonum_gcd(C_word **p, C_word n, C_word x, C_word y)
675{
676   double xub, yub, r;
677   
678   xub = C_flonum_magnitude(x);
679   yub = C_flonum_magnitude(y);
680
681   if (xub < 0.0) xub = -xub;
682   if (yub < 0.0) yub = -yub;
683   
684   while(yub != 0.0) {
685     r = fmod(xub, yub);
686     xub = yub;
687     yub = r;
688   }
689   C_return(C_flonum(p, xub));
690}
691
692
693static void
694big_plus_big(C_word c, C_word self, C_word k, C_word x, C_word y)
695{
696  bignum_type result;
697  result = bignum_add(big_of(x), big_of(y));
698  C_return_bignum(k, result);
699}
700
701/*
702 * This now makes the assumption it is never passed a bignum of LENGTH 0.
703 * This should always be valid in Chicken.
704 */
705static bignum_type
706bignum_add(bignum_type x, bignum_type y)
707{
708  return
709    (((BIGNUM_NEGATIVE_P (x))
710      ? ((BIGNUM_NEGATIVE_P (y))
711         ? (bignum_add_unsigned (x, y, 1))
712         : (bignum_subtract_unsigned (y, x)))
713      : ((BIGNUM_NEGATIVE_P (y))
714         ? (bignum_subtract_unsigned (x, y))
715         : (bignum_add_unsigned (x, y, 0)))));
716}
717
718static bignum_type
719bignum_add_unsigned(bignum_type x, bignum_type y, int negative_p)
720{
721  if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
722    {
723      bignum_type z = x;
724      x = y;
725      y = z;
726    }
727  {
728    bignum_length_type x_length = (BIGNUM_LENGTH (x));
729    bignum_type r = (bignum_allocate ((x_length + 1), negative_p));
730    bignum_digit_type sum;
731    bignum_digit_type carry = 0;
732    bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
733    bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
734    {
735      bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
736      bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
737      while (scan_y < end_y)
738        {
739          sum = ((*scan_x++) + (*scan_y++) + carry);
740          if (sum < BIGNUM_RADIX)
741            {
742              (*scan_r++) = sum;
743              carry = 0;
744            }
745          else
746            {
747              (*scan_r++) = (sum - BIGNUM_RADIX);
748              carry = 1;
749            }
750        }
751    }
752    {
753      bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
754      if (carry != 0)
755        while (scan_x < end_x)
756          {
757            sum = ((*scan_x++) + 1);
758            if (sum < BIGNUM_RADIX)
759              {
760                (*scan_r++) = sum;
761                carry = 0;
762                break;
763              }
764            else
765              (*scan_r++) = (sum - BIGNUM_RADIX);
766          }
767      while (scan_x < end_x)
768        (*scan_r++) = (*scan_x++);
769    }
770    if (carry != 0) {
771      (*scan_r) = 1;
772    } else { /* r is one word too big (to hold a possible carry), readjust */
773      BIGNUM_SET_HEADER(r, x_length, negative_p);
774      BIGNUM_REDUCE_LENGTH(r, r, x_length);
775    }
776    return (r);
777  }
778}
779
780static void
781big_minus_big(C_word c, C_word self, C_word k, C_word x, C_word y)
782{
783  bignum_type result;
784  result = bignum_subtract(big_of(x), big_of(y));
785  C_return_bignum(k, result);
786}
787
788/*
789 * This now makes the assumption it is never passed a bignum of LENGTH 0.
790 * This should always be valid in Chicken.
791 */
792static bignum_type
793bignum_subtract(bignum_type x, bignum_type y)
794{
795  return
796    (((BIGNUM_NEGATIVE_P (x))
797       ? ((BIGNUM_NEGATIVE_P (y))
798          ? (bignum_subtract_unsigned (y, x))
799          : (bignum_add_unsigned (x, y, 1)))
800       : ((BIGNUM_NEGATIVE_P (y))
801          ? (bignum_add_unsigned (x, y, 0))
802          : (bignum_subtract_unsigned (x, y)))));
803}
804
805static bignum_type
806bignum_subtract_unsigned(bignum_type x, bignum_type y)
807{
808  int negative_p;
809  switch (bignum_compare_unsigned (x, y))
810    {
811    case bignum_comparison_equal:
812      return (BIGNUM_ZERO ());
813    case bignum_comparison_less:
814      {
815        bignum_type z = x;
816        x = y;
817        y = z;
818      }
819      negative_p = 1;
820      break;
821    case bignum_comparison_greater:
822      negative_p = 0;
823      break;
824    }
825  {
826    bignum_length_type x_length = (BIGNUM_LENGTH (x));
827    bignum_type r = (bignum_allocate (x_length, negative_p));
828    bignum_digit_type difference;
829    bignum_digit_type borrow = 0;
830    bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
831    bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
832    {
833      bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
834      bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
835      while (scan_y < end_y)
836        {
837          difference = (((*scan_x++) - (*scan_y++)) - borrow);
838          if (difference < 0)
839            {
840              (*scan_r++) = (difference + BIGNUM_RADIX);
841              borrow = 1;
842            }
843          else
844            {
845              (*scan_r++) = difference;
846              borrow = 0;
847            }
848        }
849    }
850    {
851      bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
852      if (borrow != 0)
853        while (scan_x < end_x)
854          {
855            difference = ((*scan_x++) - borrow);
856            if (difference < 0)
857              (*scan_r++) = (difference + BIGNUM_RADIX);
858            else
859              {
860                (*scan_r++) = difference;
861                borrow = 0;
862                break;
863              }
864          }
865      assert(borrow == 0);
866      while (scan_x < end_x)
867        (*scan_r++) = (*scan_x++);
868    }
869    return (bignum_trim (r));
870  }
871}
872
873
874static void
875big_times_big(C_word c, C_word self, C_word k, C_word x, C_word y)
876{
877  bignum_type bigx = big_of(x), bigy = big_of(y);
878  int neg_p = ((BIGNUM_NEGATIVE_P (bigx))
879               ? (! (BIGNUM_NEGATIVE_P (bigy)))
880               : (BIGNUM_NEGATIVE_P (bigy)));
881  /* If length 1 or 0, it should be a fixnum */
882  assert(BIGNUM_LENGTH(bigx) > 1);
883  assert(BIGNUM_LENGTH(bigy) > 1);
884  C_return_bignum(k, bignum_multiply_unsigned(bigx, bigy, neg_p));
885}
886
887static void
888fix_times_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
889{
890  bignum_type bigx, bigy, result;
891  bignum_digit_type absx, absy;
892  C_word neg_p;
893
894  absx = C_unfix(x);
895  absx = absx < 0 ? -absx : absx;
896  absy = C_unfix(y);
897  absy = absy < 0 ? -absy : absy;
898  neg_p = ((x & C_INT_SIGN_BIT) ? !(y & C_INT_SIGN_BIT) : (y & C_INT_SIGN_BIT));
899
900  if (absx < BIGNUM_RADIX_ROOT) {
901     if (absx == 0 || absx == 1 || absy < BIGNUM_RADIX_ROOT) {
902       C_kontinue(k, C_fix(neg_p ? -(absx * absy) : (absx * absy)));
903     } else {
904       bigy = bignum_allocate_from_fixnum(y);
905       result = bignum_multiply_unsigned_small_factor(bigy, absx, neg_p ? 1 : 0);
906       BIGNUM_DEALLOCATE(bigy);
907       C_return_bignum(k, result);
908     }
909  } else if (absy < BIGNUM_RADIX_ROOT) {
910     if (absy == 0 || absy == 1 /*|| absx < BIGNUM_RADIX_ROOT */) {
911       C_kontinue(k, C_fix(neg_p ? -(absx * absy) : (absx * absy)));
912     } else {
913       bigx = bignum_allocate_from_fixnum(x);
914       result = bignum_multiply_unsigned_small_factor(bigx, absy, neg_p ? 1 : 0);
915       BIGNUM_DEALLOCATE(bigx);
916       C_return_bignum(k, result);
917     }
918  } else {
919    bigx = bignum_allocate_from_fixnum(x);
920    bigy = bignum_allocate_from_fixnum(y);
921    result = bignum_multiply_unsigned (bigx, bigy, neg_p ? 1 : 0);
922    BIGNUM_DEALLOCATE(bigy);
923    BIGNUM_DEALLOCATE(bigx);
924    C_return_bignum(k, result);
925  }
926}
927
928
929static void
930fix_times_big(C_word c, C_word self, C_word k, C_word x, C_word y)
931{
932  bignum_type bigx, result, bigy;
933  bignum_digit_type absx;
934  C_word neg_p;
935 
936  bigy = big_of(y);
937 
938  if (x == C_fix(0))
939    C_kontinue(k, C_fix(0));
940  else if (x == C_fix(1))
941    C_kontinue(k, y);
942  else if (x == C_fix(-1))
943    C_return_bignum(k, bignum_new_sign(bigy, !(BIGNUM_NEGATIVE_P(bigy))));
944
945  absx = C_unfix(x);
946  absx = absx < 0 ? -absx : absx; 
947  neg_p = ((x & C_INT_SIGN_BIT)
948           ? !(BIGNUM_NEGATIVE_P(bigy)) : (BIGNUM_NEGATIVE_P(bigy)));
949 
950  if (absx < BIGNUM_RADIX_ROOT) {
951     C_return_bignum(k, bignum_multiply_unsigned_small_factor(bigy, absx, neg_p));
952  } else {
953    bigx = bignum_allocate_from_fixnum(x);
954    result = bignum_multiply_unsigned (bigx, bigy, neg_p);
955    BIGNUM_DEALLOCATE(bigx);
956    C_return_bignum(k, result);
957  }
958}
959
960/* Multiplication
961   Maximum value for product_low or product_high:
962        ((R * R) + (R * (R - 2)) + (R - 1))
963   Maximum value for carry: ((R * (R - 1)) + (R - 1))
964        where R == BIGNUM_RADIX_ROOT */
965
966static bignum_type
967bignum_multiply_unsigned(bignum_type x, bignum_type y, int negative_p)
968{
969  if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
970    {
971      bignum_type z = x;
972      x = y;
973      y = z;
974    }
975  {
976    bignum_digit_type carry;
977    bignum_digit_type y_digit_low;
978    bignum_digit_type y_digit_high;
979    bignum_digit_type x_digit_low;
980    bignum_digit_type x_digit_high;
981    bignum_digit_type product_low;
982    bignum_digit_type * scan_r;
983    bignum_digit_type * scan_y;
984    bignum_length_type x_length = (BIGNUM_LENGTH (x));
985    bignum_length_type y_length = (BIGNUM_LENGTH (y));
986    bignum_type r = (bignum_allocate_zeroed ((x_length + y_length), negative_p));
987    bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
988    bignum_digit_type * end_x = (scan_x + x_length);
989    bignum_digit_type * start_y = (BIGNUM_START_PTR (y));
990    bignum_digit_type * end_y = (start_y + y_length);
991    bignum_digit_type * start_r = (BIGNUM_START_PTR (r));
992#define x_digit x_digit_high
993#define y_digit y_digit_high
994#define product_high carry
995    while (scan_x < end_x)
996      {
997        x_digit = (*scan_x++);
998        x_digit_low = (HD_LOW (x_digit));
999        x_digit_high = (HD_HIGH (x_digit));
1000        carry = 0;
1001        scan_y = start_y;
1002        scan_r = (start_r++);
1003        while (scan_y < end_y)
1004          {
1005            y_digit = (*scan_y++);
1006            y_digit_low = (HD_LOW (y_digit));
1007            y_digit_high = (HD_HIGH (y_digit));
1008            product_low =
1009              ((*scan_r) +
1010               (x_digit_low * y_digit_low) +
1011               (HD_LOW (carry)));
1012            product_high =
1013              ((x_digit_high * y_digit_low) +
1014               (x_digit_low * y_digit_high) +
1015               (HD_HIGH (product_low)) +
1016               (HD_HIGH (carry)));
1017            (*scan_r++) =
1018              (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
1019            carry =
1020              ((x_digit_high * y_digit_high) +
1021               (HD_HIGH (product_high)));
1022          }
1023        (*scan_r) += carry;
1024      }
1025    return (bignum_trim (r));
1026#undef x_digit
1027#undef y_digit
1028#undef product_high
1029  }
1030}
1031
1032static bignum_type
1033bignum_multiply_unsigned_small_factor(bignum_type x, bignum_digit_type y,
1034                                      int negative_p)
1035{
1036  bignum_length_type length_x = (BIGNUM_LENGTH (x));
1037  bignum_type p = (bignum_allocate ((length_x + 1), negative_p));
1038  bignum_destructive_copy (x, p);
1039  (BIGNUM_REF (p, length_x)) = 0;
1040  bignum_destructive_scale_up (p, y);
1041  return (bignum_trim (p));
1042}
1043
1044static void
1045bignum_destructive_scale_up(bignum_type bignum, bignum_digit_type factor)
1046{
1047  bignum_digit_type carry = 0;
1048  bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
1049  bignum_digit_type two_digits;
1050  bignum_digit_type product_low;
1051#define product_high carry
1052  bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum)));
1053  assert((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
1054  while (scan < end)
1055    {
1056      two_digits = (*scan);
1057      product_low = ((factor * (HD_LOW (two_digits))) + (HD_LOW (carry)));
1058      product_high =
1059        ((factor * (HD_HIGH (two_digits))) +
1060         (HD_HIGH (product_low)) +
1061         (HD_HIGH (carry)));
1062      (*scan++) = (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
1063      carry = (HD_HIGH (product_high));
1064    }
1065  /* A carry here would be an overflow, i.e. it would not fit.
1066     Hopefully the callers allocate enough space that this will
1067     never happen.
1068   */
1069  assert(carry == 0);
1070  return;
1071#undef product_high
1072}
1073
1074static void
1075bignum_destructive_add(bignum_type bignum, bignum_digit_type n)
1076{
1077  bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
1078  bignum_digit_type digit;
1079  digit = ((*scan) + n);
1080  if (digit < BIGNUM_RADIX)
1081    {
1082      (*scan) = digit;
1083      return;
1084    }
1085  (*scan++) = (digit - BIGNUM_RADIX);
1086  while (1)
1087    {
1088      digit = ((*scan) + 1);
1089      if (digit < BIGNUM_RADIX)
1090        {
1091          (*scan) = digit;
1092          return;
1093        }
1094      (*scan++) = (digit - BIGNUM_RADIX);
1095    }
1096}
1097
1098/* Division */
1099
1100/* For help understanding this algorithm, see:
1101   Knuth, Donald E., "The Art of Computer Programming",
1102   volume 2, "Seminumerical Algorithms"
1103   section 4.3.1, "Multiple-Precision Arithmetic". */
1104
1105static void
1106bignum_divide_unsigned_large_denominator(bignum_type numerator,
1107                                         bignum_type denominator,
1108                                         bignum_type * quotient,
1109                                         bignum_type * remainder,
1110                                         int q_negative_p,
1111                                         int r_negative_p)
1112{
1113  bignum_length_type length_n = ((BIGNUM_LENGTH (numerator)) + 1);
1114  bignum_length_type length_d = (BIGNUM_LENGTH (denominator));
1115  bignum_type q =
1116    ((quotient != ((bignum_type *) 0))
1117     ? (bignum_allocate ((length_n - length_d), q_negative_p))
1118     : BIGNUM_OUT_OF_BAND);
1119  bignum_type u = (bignum_allocate (length_n, r_negative_p));
1120  int shift = 0;
1121  assert(length_d > 1);
1122  {
1123    bignum_digit_type v1 = (BIGNUM_REF ((denominator), (length_d - 1)));
1124    while (v1 < (BIGNUM_RADIX / 2))
1125      {
1126        v1 <<= 1;
1127        shift += 1;
1128      }
1129  }
1130  if (shift == 0)
1131    {
1132      bignum_destructive_copy (numerator, u);
1133      (BIGNUM_REF (u, (length_n - 1))) = 0;
1134      bignum_divide_unsigned_normalized (u, denominator, q);
1135      if (remainder != ((bignum_type *) 0))
1136        (*remainder) = (bignum_trim (u));
1137      else
1138        BIGNUM_DEALLOCATE (u);
1139    }
1140  else
1141    {
1142      bignum_type v = (bignum_allocate (length_d, 0));
1143      bignum_destructive_normalization (numerator, u, shift);
1144      bignum_destructive_normalization (denominator, v, shift);
1145      bignum_divide_unsigned_normalized (u, v, q);
1146      BIGNUM_DEALLOCATE (v);
1147      if (remainder != ((bignum_type *) 0))
1148        (*remainder) = bignum_destructive_unnormalization (u, shift);
1149      else
1150        BIGNUM_DEALLOCATE(u);
1151    }
1152  if (quotient != ((bignum_type *) 0))
1153    (*quotient) = (bignum_trim (q));
1154  return;
1155}
1156
1157static void
1158bignum_divide_unsigned_normalized(bignum_type u, bignum_type v, bignum_type q)
1159{
1160  bignum_length_type u_length = (BIGNUM_LENGTH (u));
1161  bignum_length_type v_length = (BIGNUM_LENGTH (v));
1162  bignum_digit_type * u_start = (BIGNUM_START_PTR (u));
1163  bignum_digit_type * u_scan = (u_start + u_length);
1164  bignum_digit_type * u_scan_limit = (u_start + v_length);
1165  bignum_digit_type * u_scan_start = (u_scan - v_length);
1166  bignum_digit_type * v_start = (BIGNUM_START_PTR (v));
1167  bignum_digit_type * v_end = (v_start + v_length);
1168  bignum_digit_type * q_scan;
1169  bignum_digit_type v1 = (v_end[-1]);
1170  bignum_digit_type v2 = (v_end[-2]);
1171  bignum_digit_type ph; /* high half of double-digit product */
1172  bignum_digit_type pl; /* low half of double-digit product */
1173  bignum_digit_type guess;
1174  bignum_digit_type gh; /* high half-digit of guess */
1175  bignum_digit_type ch; /* high half of double-digit comparand */
1176  bignum_digit_type v2l = (HD_LOW (v2));
1177  bignum_digit_type v2h = (HD_HIGH (v2));
1178  bignum_digit_type cl; /* low half of double-digit comparand */
1179#define gl ph                   /* low half-digit of guess */
1180#define uj pl
1181#define qj ph
1182  bignum_digit_type gm;         /* memory loc for reference parameter */
1183  if (q != BIGNUM_OUT_OF_BAND)
1184    q_scan = ((BIGNUM_START_PTR (q)) + (BIGNUM_LENGTH (q)));
1185  while (u_scan_limit < u_scan)
1186    {
1187      uj = (*--u_scan);
1188      if (uj != v1)
1189        {
1190          /* comparand =
1191             (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
1192             guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
1193          cl = (u_scan[-2]);
1194          ch = (bignum_digit_divide (uj, (u_scan[-1]), v1, (&gm)));
1195          guess = gm;
1196        }
1197      else
1198        {
1199          cl = (u_scan[-2]);
1200          ch = ((u_scan[-1]) + v1);
1201          guess = (BIGNUM_RADIX - 1);
1202        }
1203      while (1)
1204        {
1205          /* product = (guess * v2); */
1206          gl = (HD_LOW (guess));
1207          gh = (HD_HIGH (guess));
1208          pl = (v2l * gl);
1209          ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH (pl)));
1210          pl = (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl))));
1211          ph = ((v2h * gh) + (HD_HIGH (ph)));
1212          /* if (comparand >= product) */
1213          if ((ch > ph) || ((ch == ph) && (cl >= pl)))
1214            break;
1215          guess -= 1;
1216          /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
1217          ch += v1;
1218          /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
1219          if (ch >= BIGNUM_RADIX)
1220            break;
1221        }
1222      qj = (bignum_divide_subtract (v_start, v_end, guess, (--u_scan_start)));
1223      if (q != BIGNUM_OUT_OF_BAND)
1224        (*--q_scan) = qj;
1225    }
1226  return;
1227#undef gl
1228#undef uj
1229#undef qj
1230}
1231
1232static bignum_digit_type
1233bignum_divide_subtract(bignum_digit_type * v_start,
1234                       bignum_digit_type * v_end,
1235                       bignum_digit_type guess,
1236                       bignum_digit_type * u_start)
1237{
1238  bignum_digit_type * v_scan = v_start;
1239  bignum_digit_type * u_scan = u_start;
1240  bignum_digit_type carry = 0;
1241  if (guess == 0) return (0);
1242  {
1243    bignum_digit_type gl = (HD_LOW (guess));
1244    bignum_digit_type gh = (HD_HIGH (guess));
1245    bignum_digit_type v;
1246    bignum_digit_type pl;
1247    bignum_digit_type vl;
1248#define vh v
1249#define ph carry
1250#define diff pl
1251    while (v_scan < v_end)
1252      {
1253        v = (*v_scan++);
1254        vl = (HD_LOW (v));
1255        vh = (HD_HIGH (v));
1256        pl = ((vl * gl) + (HD_LOW (carry)));
1257        ph = ((vl * gh) + (vh * gl) + (HD_HIGH (pl)) + (HD_HIGH (carry)));
1258        diff = ((*u_scan) - (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl)))));
1259        if (diff < 0)
1260          {
1261            (*u_scan++) = (diff + BIGNUM_RADIX);
1262            carry = ((vh * gh) + (HD_HIGH (ph)) + 1);
1263          }
1264        else
1265          {
1266            (*u_scan++) = diff;
1267            carry = ((vh * gh) + (HD_HIGH (ph)));
1268          }
1269      }
1270    if (carry == 0)
1271      return (guess);
1272    diff = ((*u_scan) - carry);
1273    if (diff < 0)
1274      (*u_scan) = (diff + BIGNUM_RADIX);
1275    else
1276      {
1277        (*u_scan) = diff;
1278        return (guess);
1279      }
1280#undef vh
1281#undef ph
1282#undef diff
1283  }
1284  /* Subtraction generated carry, implying guess is one too large.
1285     Add v back in to bring it back down. */
1286  v_scan = v_start;
1287  u_scan = u_start;
1288  carry = 0;
1289  while (v_scan < v_end)
1290    {
1291      bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
1292      if (sum < BIGNUM_RADIX)
1293        {
1294          (*u_scan++) = sum;
1295          carry = 0;
1296        }
1297      else
1298        {
1299          (*u_scan++) = (sum - BIGNUM_RADIX);
1300          carry = 1;
1301        }
1302    }
1303  if (carry == 1)
1304    {
1305      bignum_digit_type sum = ((*u_scan) + carry);
1306      (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
1307    }
1308  return (guess - 1);
1309}
1310
1311static void
1312bignum_divide_unsigned_medium_denominator(bignum_type numerator,
1313                                          bignum_digit_type denominator,
1314                                          bignum_type * quotient,
1315                                          bignum_type * remainder,
1316                                          int q_negative_p,
1317                                          int r_negative_p)
1318{
1319  bignum_length_type length_n = (BIGNUM_LENGTH (numerator));
1320  bignum_length_type length_q;
1321  bignum_type q;
1322  int shift = 0;
1323  /* Because `bignum_digit_divide' requires a normalized denominator. */
1324  while (denominator < (BIGNUM_RADIX / 2))
1325    {
1326      denominator <<= 1;
1327      shift += 1;
1328    }
1329  if (shift == 0)
1330    {
1331      length_q = length_n;
1332      q = (bignum_allocate (length_q, q_negative_p));
1333      bignum_destructive_copy (numerator, q);
1334    }
1335  else
1336    {
1337      length_q = (length_n + 1);
1338      q = (bignum_allocate (length_q, q_negative_p));
1339      bignum_destructive_normalization (numerator, q, shift);
1340    }
1341  {
1342    bignum_digit_type r = 0;
1343    bignum_digit_type * start = (BIGNUM_START_PTR (q));
1344    bignum_digit_type * scan = (start + length_q);
1345    bignum_digit_type qj;
1346    if (quotient != ((bignum_type *) 0))
1347      {
1348        while (start < scan)
1349          {
1350            r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
1351            (*scan) = qj;
1352          }
1353        (*quotient) = (bignum_trim (q));
1354      }
1355    else
1356      {
1357        while (start < scan)
1358          r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
1359        BIGNUM_DEALLOCATE (q);
1360      }
1361    if (remainder != ((bignum_type *) 0))
1362      {
1363        if (shift != 0)
1364          r >>= shift;
1365        (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
1366      }
1367  }
1368  return;
1369}
1370
1371static void
1372bignum_destructive_normalization(bignum_type source, bignum_type target,
1373                                 int shift_left)
1374{
1375  bignum_digit_type digit;
1376  bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
1377  bignum_digit_type carry = 0;
1378  bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
1379  bignum_digit_type * end_source = (scan_source + (BIGNUM_LENGTH (source)));
1380  bignum_digit_type * end_target = (scan_target + (BIGNUM_LENGTH (target)));
1381  int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
1382  bignum_digit_type mask = ((1L << shift_right) - 1);
1383  while (scan_source < end_source)
1384    {
1385      digit = (*scan_source++);
1386      (*scan_target++) = (((digit & mask) << shift_left) | carry);
1387      carry = (digit >> shift_right);
1388    }
1389  if (scan_target < end_target)
1390    (*scan_target) = carry;
1391  else
1392    assert(carry == 0);
1393  return;
1394}
1395
1396/* This will also trim the number if necessary */
1397static bignum_type
1398bignum_destructive_unnormalization(bignum_type bignum, int shift_right)
1399{
1400  bignum_length_type length = BIGNUM_LENGTH(bignum);
1401  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
1402  bignum_digit_type * scan = (start + length);
1403  bignum_digit_type digit;
1404  bignum_digit_type carry = 0;
1405  int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
1406  bignum_digit_type mask = ((1L << shift_right) - 1);
1407 
1408  while (!(*--scan)) {
1409    if (start == scan) { /* Don't bother. */
1410      BIGNUM_SET_HEADER (bignum, 0, 0);
1411      BIGNUM_REDUCE_LENGTH (bignum, bignum, 0);
1412      return bignum;
1413    }
1414    --length;
1415  }
1416
1417  digit = (*scan);
1418  (*scan) = (digit >> shift_right);
1419  length -= (*scan == 0); /* Add 1 or 0 */
1420  carry = ((digit & mask) << shift_left);
1421 
1422  while (start < scan)
1423    {
1424      digit = (*--scan);
1425      (*scan) = ((digit >> shift_right) | carry);
1426      carry = ((digit & mask) << shift_left);
1427    }
1428  assert(carry == 0);
1429  assert(BIGNUM_LENGTH(bignum) != length);
1430  assert(length != 1 || BIGNUM_REF(bignum, 0) != 0);
1431  BIGNUM_SET_HEADER
1432    (bignum, length, (BIGNUM_NEGATIVE_P (bignum)));
1433  BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
1434
1435  return bignum;
1436}
1437
1438/* This is a reduced version of the division algorithm, applied to the
1439   case of dividing two bignum digits by one bignum digit.  It is
1440   assumed that the numerator, denominator are normalized. */
1441
1442#define BDD_STEP(qn, j)                                                 \
1443{                                                                       \
1444  uj = (u[j]);                                                          \
1445  if (uj != v1)                                                         \
1446    {                                                                   \
1447      uj_uj1 = (HD_CONS (uj, (u[j + 1])));                              \
1448      guess = (uj_uj1 / v1);                                            \
1449      comparand = (HD_CONS ((uj_uj1 % v1), (u[j + 2])));                \
1450    }                                                                   \
1451  else                                                                  \
1452    {                                                                   \
1453      guess = (BIGNUM_RADIX_ROOT - 1);                                  \
1454      comparand = (HD_CONS (((u[j + 1]) + v1), (u[j + 2])));            \
1455    }                                                                   \
1456  while ((guess * v2) > comparand)                                      \
1457    {                                                                   \
1458      guess -= 1;                                                       \
1459      comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH);                    \
1460      if (comparand >= BIGNUM_RADIX)                                    \
1461        break;                                                          \
1462    }                                                                   \
1463  qn = (bignum_digit_divide_subtract (v1, v2, guess, (&u[j])));         \
1464}
1465
1466static bignum_digit_type
1467bignum_digit_divide(bignum_digit_type uh, bignum_digit_type ul,
1468                    bignum_digit_type v,
1469                    bignum_digit_type * q) /* return value */
1470{
1471  bignum_digit_type guess;
1472  bignum_digit_type comparand;
1473  bignum_digit_type v1;
1474  bignum_digit_type v2;
1475  bignum_digit_type uj;
1476  bignum_digit_type uj_uj1;
1477  bignum_digit_type q1;
1478  bignum_digit_type q2;
1479  bignum_digit_type u [4];
1480  if (uh == 0)
1481    {
1482      if (ul < v)
1483        {
1484          (*q) = 0;
1485          return (ul);
1486        }
1487      else if (ul == v)
1488        {
1489          (*q) = 1;
1490          return (0);
1491        }
1492    }
1493  (u[0]) = (HD_HIGH (uh));
1494  (u[1]) = (HD_LOW (uh));
1495  (u[2]) = (HD_HIGH (ul));
1496  (u[3]) = (HD_LOW (ul));
1497  v1 = (HD_HIGH (v));
1498  v2 = (HD_LOW (v));
1499  BDD_STEP (q1, 0);
1500  BDD_STEP (q2, 1);
1501  (*q) = (HD_CONS (q1, q2));
1502  return (HD_CONS ((u[2]), (u[3])));
1503}
1504
1505#undef BDD_STEP
1506
1507#define BDDS_MULSUB(vn, un, carry_in)                                   \
1508{                                                                       \
1509  product = ((vn * guess) + carry_in);                                  \
1510  diff = (un - (HD_LOW (product)));                                     \
1511  if (diff < 0)                                                         \
1512    {                                                                   \
1513      un = (diff + BIGNUM_RADIX_ROOT);                                  \
1514      carry = ((HD_HIGH (product)) + 1);                                \
1515    }                                                                   \
1516  else                                                                  \
1517    {                                                                   \
1518      un = diff;                                                        \
1519      carry = (HD_HIGH (product));                                      \
1520    }                                                                   \
1521}
1522
1523#define BDDS_ADD(vn, un, carry_in)                                      \
1524{                                                                       \
1525  sum = (vn + un + carry_in);                                           \
1526  if (sum < BIGNUM_RADIX_ROOT)                                          \
1527    {                                                                   \
1528      un = sum;                                                         \
1529      carry = 0;                                                        \
1530    }                                                                   \
1531  else                                                                  \
1532    {                                                                   \
1533      un = (sum - BIGNUM_RADIX_ROOT);                                   \
1534      carry = 1;                                                        \
1535    }                                                                   \
1536}
1537
1538static bignum_digit_type
1539bignum_digit_divide_subtract(bignum_digit_type v1, bignum_digit_type v2,
1540                             bignum_digit_type guess, bignum_digit_type * u)
1541{
1542  {
1543    bignum_digit_type product;
1544    bignum_digit_type diff;
1545    bignum_digit_type carry;
1546    BDDS_MULSUB (v2, (u[2]), 0);
1547    BDDS_MULSUB (v1, (u[1]), carry);
1548    if (carry == 0)
1549      return (guess);
1550    diff = ((u[0]) - carry);
1551    if (diff < 0)
1552      (u[0]) = (diff + BIGNUM_RADIX);
1553    else
1554      {
1555        (u[0]) = diff;
1556        return (guess);
1557      }
1558  }
1559  {
1560    bignum_digit_type sum;
1561    bignum_digit_type carry;
1562    BDDS_ADD(v2, (u[2]), 0);
1563    BDDS_ADD(v1, (u[1]), carry);
1564    if (carry == 1)
1565      (u[0]) += 1;
1566  }
1567  return (guess - 1);
1568}
1569
1570#undef BDDS_MULSUB
1571#undef BDDS_ADD
1572
1573static void
1574bignum_divide_unsigned_small_denominator(bignum_type numerator,
1575                                         bignum_digit_type denominator,
1576                                         bignum_type * quotient,
1577                                         bignum_type * remainder,
1578                                         int q_negative_p,  int r_negative_p)
1579{
1580  bignum_type q = (bignum_new_sign (numerator, q_negative_p));
1581  bignum_digit_type r = (bignum_destructive_scale_down (q, denominator));
1582  (*quotient) = (bignum_trim (q));
1583  if (remainder != ((bignum_type *) 0))
1584    (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
1585  return;
1586}
1587
1588/* Given (denominator > 1), it is fairly easy to show that
1589   (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
1590   that all digits are < BIGNUM_RADIX. */
1591
1592static bignum_digit_type
1593bignum_destructive_scale_down(bignum_type bignum, bignum_digit_type denominator)
1594{
1595  bignum_digit_type numerator;
1596  bignum_digit_type remainder = 0;
1597  bignum_digit_type two_digits;
1598#define quotient_high remainder
1599  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
1600  bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
1601  assert((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
1602  while (start < scan)
1603    {
1604      two_digits = (*--scan);
1605      numerator = (HD_CONS (remainder, (HD_HIGH (two_digits))));
1606      quotient_high = (numerator / denominator);
1607      numerator = (HD_CONS ((numerator % denominator), (HD_LOW (two_digits))));
1608      (*scan) = (HD_CONS (quotient_high, (numerator / denominator)));
1609      remainder = (numerator % denominator);
1610    }
1611  return (remainder);
1612#undef quotient_high
1613}
1614
1615static bignum_type
1616bignum_remainder_unsigned_small_denominator(bignum_type n, bignum_digit_type d,
1617                                            int negative_p)
1618{
1619  bignum_digit_type two_digits;
1620  bignum_digit_type * start = (BIGNUM_START_PTR (n));
1621  bignum_digit_type * scan = (start + (BIGNUM_LENGTH (n)));
1622  bignum_digit_type r = 0;
1623  assert((d > 1) && (d < BIGNUM_RADIX_ROOT));
1624  while (start < scan)
1625    {
1626      two_digits = (*--scan);
1627      r =
1628        ((HD_CONS (((HD_CONS (r, (HD_HIGH (two_digits)))) % d),
1629                   (HD_LOW (two_digits))))
1630         % d);
1631    }
1632  return (bignum_digit_to_bignum (r, negative_p));
1633}
1634
1635static void
1636fix_neg(C_word c, C_word self, C_word k, C_word x)
1637{
1638  x = C_unfix(x);
1639  /* This code "knows" that bignums have 2 "reserved" bits, like fixnums */
1640  if (x != C_MOST_NEGATIVE_FIXNUM) { /* C_fitsinfixnump(x) */
1641    C_kontinue(k, C_fix(-x));
1642  } else {
1643    bignum_digit_type *ref;
1644    bignum_type big;
1645    big = bignum_allocate(2, 0);
1646    ref = BIGNUM_START_PTR(big);
1647    *ref++ = 0;
1648    *ref = 1;
1649    C_return_bignum(k, big);
1650  }
1651}
1652
1653static void
1654big_neg(C_word c, C_word self, C_word k, C_word x)
1655{
1656  bignum_type big = big_of(x);
1657  C_return_bignum(k, bignum_new_sign(big, !(BIGNUM_NEGATIVE_P (big))));
1658}
1659
1660static C_word
1661big_comp_big(C_word x, C_word y)
1662{
1663  bignum_type bigx = big_of(x), bigy = big_of(y);
1664  return C_fix((BIGNUM_NEGATIVE_P (bigx))
1665               ? ((BIGNUM_NEGATIVE_P (bigy))
1666                  ? (bignum_compare_unsigned (bigy, bigx))
1667                  : (bignum_comparison_less))
1668               : ((BIGNUM_NEGATIVE_P (bigy))
1669                  ? (bignum_comparison_greater)
1670                  : (bignum_compare_unsigned (bigx, bigy))));
1671}
1672
1673static enum bignum_comparison
1674bignum_compare_unsigned(bignum_type x, bignum_type y)
1675{
1676  bignum_length_type x_length;
1677  bignum_length_type y_length;
1678  if (x == y) /* Objects are the same? */
1679    return (bignum_comparison_equal);
1680
1681  x_length = (BIGNUM_LENGTH (x));
1682  y_length = (BIGNUM_LENGTH (y));
1683  if (x_length < y_length)
1684    return (bignum_comparison_less);
1685  if (x_length > y_length)
1686    return (bignum_comparison_greater);
1687  {
1688    bignum_digit_type * start_x = (BIGNUM_START_PTR (x));
1689    bignum_digit_type * scan_x = (start_x + x_length);
1690    bignum_digit_type * scan_y = ((BIGNUM_START_PTR (y)) + y_length);
1691    while (start_x < scan_x)
1692      {
1693        bignum_digit_type digit_x = (*--scan_x);
1694        bignum_digit_type digit_y = (*--scan_y);
1695        if (digit_x < digit_y)
1696          return (bignum_comparison_less);
1697        if (digit_x > digit_y)
1698          return (bignum_comparison_greater);
1699      }
1700  }
1701  return (bignum_comparison_equal);
1702}
1703
1704static void
1705big_abs(C_word c, C_word self, C_word k, C_word big)
1706{
1707  if (!BIGNUM_NEGATIVE_P(big_of(big)))
1708    C_kontinue(k, big);
1709  else
1710    C_return_bignum(k, bignum_new_sign(big_of(big), 0));
1711}
1712
1713static void
1714big_quotient_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
1715{
1716  bignum_type bigx = big_of(x);
1717  y = C_unfix(y);
1718
1719  if (y == 1)
1720    C_kontinue(k, x);
1721  else if (y == -1)
1722    C_return_bignum(k, bignum_new_sign(bigx, !(BIGNUM_NEGATIVE_P(bigx))));
1723
1724  /* Too bad, we really need to do some work... */
1725  {
1726    int neg_p = (y < 0) ? !(BIGNUM_NEGATIVE_P(bigx)) : BIGNUM_NEGATIVE_P(bigx);
1727    bignum_digit_type abs_y = (y < 0) ? -y : y;
1728    bignum_type quotient;
1729   
1730    if (y == C_MOST_NEGATIVE_FIXNUM) {
1731      if (!BIGNUM_NEGATIVE_P(bigx) && BIGNUM_LENGTH(bigx) == 1
1732          && BIGNUM_REF(bigx, 1) == 1 && BIGNUM_REF(bigx, 0) == 0) {
1733        /*
1734         * Very very special case:
1735         * quotient(MOST_NEGATIVE_FIXNUM, -(MOST_NEGATIVE_FIXNUM)) => -1
1736         */
1737        C_kontinue(k, C_fix(-1));
1738      } else {
1739        /* This is the only case we need to go allocate a bignum for */
1740        bignum_type bigy =
1741          bignum_allocate_from_fixnum(C_fix(C_MOST_NEGATIVE_FIXNUM));
1742
1743        bignum_divide_unsigned_large_denominator
1744          (bigx, bigy, (&quotient), ((bignum_type *) 0), neg_p, 0);
1745        BIGNUM_DEALLOCATE(bigy);
1746        C_return_bignum(k, quotient);
1747      }
1748    } else if (abs_y < BIGNUM_RADIX_ROOT) {
1749      bignum_divide_unsigned_small_denominator
1750        (bigx, abs_y, (&quotient), ((bignum_type *) 0), neg_p, 0);
1751      C_return_bignum(k, quotient);
1752    } else {
1753      bignum_divide_unsigned_medium_denominator
1754        (bigx, abs_y, (&quotient), ((bignum_type *) 0), neg_p, 0);
1755      C_return_bignum(k, quotient);
1756    }
1757  }
1758}
1759
1760static void
1761big_quotient_big(C_word c, C_word self, C_word k, C_word x, C_word y)
1762{
1763  bignum_type numerator = big_of(x);
1764  bignum_type denominator = big_of(y);
1765  int neg_p = ((BIGNUM_NEGATIVE_P (denominator))
1766               ? (! (BIGNUM_NEGATIVE_P (numerator)))
1767               : (BIGNUM_NEGATIVE_P (numerator)));
1768
1769  switch (bignum_compare_unsigned (numerator, denominator))
1770    {
1771    case bignum_comparison_equal:
1772      C_kontinue(k, neg_p ? C_fix(-1) : C_fix(1));
1773    case bignum_comparison_less:
1774      C_kontinue(k, C_fix(0));
1775    case bignum_comparison_greater:
1776    default:                                  /* to appease gcc -Wall */
1777      {
1778        bignum_type quotient;
1779        assert(BIGNUM_LENGTH(denominator) > 1);
1780        bignum_divide_unsigned_large_denominator
1781          (numerator, denominator, (&quotient), ((bignum_type *) 0), neg_p, 0);
1782        C_return_bignum(k, quotient);
1783      }
1784    }
1785}
1786
1787static void
1788big_remainder_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
1789{
1790   bignum_type remainder, bigy, bigx;
1791   int x_neg, y_neg;
1792   C_word abs_y;
1793
1794   switch (y) {
1795   /* case 0:  SHOULD NOT HAPPEN (checked in Scheme)
1796     C_kontinue(k, C_SCHEME_FALSE); */
1797   case C_fix(1):
1798   case C_fix(-1):
1799     C_kontinue(k, C_fix(0));
1800   case C_fix(C_MOST_NEGATIVE_FIXNUM):
1801     bigx = big_of(x);
1802     if (BIGNUM_LENGTH(bigx) == 2 && BIGNUM_REF(bigx, 0) == 0 && BIGNUM_REF(bigx, 0))
1803       C_kontinue(k, C_fix(0));
1804     /* Don't handle 0 <= length(bigx) <= 1 since then it should be a fixnum */
1805     assert(BIGNUM_LENGTH(bigx) >= 2);
1806     
1807     bigy = bignum_allocate_from_fixnum(y);
1808     bignum_divide_unsigned_large_denominator
1809      (bigx, bigy, (bignum_type *)0, &remainder, 0, BIGNUM_NEGATIVE_P(bigx));
1810     BIGNUM_DEALLOCATE(bigy);
1811     
1812     C_return_bignum(k, remainder);
1813   default:
1814     bigx = big_of(x);
1815     y = C_unfix(y);
1816     y_neg = (y < 0);
1817     abs_y = y_neg ? -y : y;
1818     x_neg = BIGNUM_NEGATIVE_P(bigx);
1819     
1820     if (abs_y < BIGNUM_RADIX_ROOT)
1821       remainder =
1822         bignum_remainder_unsigned_small_denominator(bigx, abs_y, x_neg);
1823     else
1824       bignum_divide_unsigned_medium_denominator(bigx, abs_y,
1825                                                 (bignum_type *) 0, &remainder,
1826                                                 x_neg, x_neg);
1827     C_return_bignum(k, remainder);
1828   }
1829}
1830
1831static void
1832big_remainder_big(C_word c, C_word self, C_word k, C_word x, C_word y)
1833{
1834  bignum_type numerator = big_of(x), denominator = big_of(y);
1835 
1836  switch (bignum_compare_unsigned (numerator, denominator))
1837    {
1838      case bignum_comparison_equal:
1839        C_kontinue(k, C_fix(0));
1840      case bignum_comparison_less:
1841        C_kontinue(k, x);
1842      case bignum_comparison_greater:
1843      default:                                  /* to appease gcc -Wall */
1844        {
1845          bignum_type remainder;
1846          bignum_divide_unsigned_large_denominator
1847            (numerator, denominator,
1848             ((bignum_type *) 0), (&remainder),
1849             0, BIGNUM_NEGATIVE_P(numerator));
1850          C_return_bignum(k, remainder);
1851        }
1852    }
1853}
1854
1855static void
1856flo_to_int(C_word c, C_word self, C_word k, C_word x)
1857{
1858  C_return_bignum(k, double_to_bignum(C_flonum_magnitude(x)));
1859}
1860
1861#define DTB_WRITE_DIGIT(factor)                                         \
1862{                                                                       \
1863  significand *= (factor);                                              \
1864  digit = ((bignum_digit_type) significand);                            \
1865  (*--scan) = digit;                                                    \
1866  significand -= ((double) digit);                                      \
1867}
1868
1869static bignum_type
1870double_to_bignum(double x)
1871{
1872  int exponent;
1873  double significand = (frexp (x, (&exponent)));
1874  if (exponent <= 0) return (BIGNUM_ZERO ());
1875  if (exponent == 1) return (BIGNUM_ONE (x < 0));
1876  if (significand < 0) significand = (-significand);
1877  {
1878    bignum_length_type length = (BIGNUM_BITS_TO_DIGITS (exponent));
1879    bignum_type result = (bignum_allocate (length, (x < 0)));
1880    bignum_digit_type * start = (BIGNUM_START_PTR (result));
1881    bignum_digit_type * scan = (start + length);
1882    bignum_digit_type digit;
1883    int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
1884    if (odd_bits > 0)
1885      DTB_WRITE_DIGIT (1L << odd_bits);
1886    while (start < scan)
1887      {
1888        if (significand == 0)
1889          {
1890            while (start < scan)
1891              (*--scan) = 0;
1892            break;
1893          }
1894        DTB_WRITE_DIGIT (BIGNUM_RADIX);
1895      }
1896    return (result);
1897  }
1898}
1899
1900#undef DTB_WRITE_DIGIT
1901
1902static void 
1903int_bitwise_int(C_word c, C_word self, C_word k, C_word op, C_word x, C_word y)
1904{
1905  bignum_type bigx, bigy, result;
1906 
1907  if((x & C_FIXNUM_BIT) != 0)
1908    bigx = bignum_allocate_from_fixnum(x);
1909  else
1910    bigx = big_of(x);
1911
1912  if((y & C_FIXNUM_BIT) != 0)
1913    bigy = bignum_allocate_from_fixnum(y);
1914  else
1915    bigy = big_of(y);
1916
1917  result = ((BIGNUM_NEGATIVE_P (bigx))
1918             ? (BIGNUM_NEGATIVE_P (bigy))
1919               ? bignum_negneg_bitwise_op(C_unfix(op), bigx, bigy)
1920               : bignum_posneg_bitwise_op(C_unfix(op), bigy, bigx)
1921             : (BIGNUM_NEGATIVE_P (bigy))
1922               ? bignum_posneg_bitwise_op(C_unfix(op), bigx, bigy)
1923               : bignum_pospos_bitwise_op(C_unfix(op), bigx, bigy));
1924       
1925  if((x & C_FIXNUM_BIT) != 0)
1926    BIGNUM_DEALLOCATE(bigx);
1927  if((y & C_FIXNUM_BIT) != 0)
1928    BIGNUM_DEALLOCATE(bigy);
1929
1930  C_return_bignum(k, result);
1931}
1932
1933static void 
1934int_not(C_word c, C_word self, C_word k, C_word x)
1935{
1936  bignum_type bigx, result;
1937 
1938  if((x & C_FIXNUM_BIT) != 0)
1939    bigx = bignum_allocate_from_fixnum(x);
1940  else
1941    bigx = big_of(x);
1942
1943  result = bignum_bitwise_not(bigx);
1944
1945  if((x & C_FIXNUM_BIT) != 0)
1946    BIGNUM_DEALLOCATE(bigx);
1947
1948  C_return_bignum(k, result);
1949}
1950
1951bignum_type
1952bignum_bitwise_not(bignum_type x)
1953{
1954  static C_word invbits[] = { BIGNUM_RADIX | 1, 1 }; /* bignum representing -1 */
1955  return bignum_subtract((bignum_type)invbits, x);
1956}
1957
1958/*
1959 * From Hacker's Delight by Frank Warren
1960 * based on a modified nlz() from section 5-3 (fig. 5-7)
1961 */
1962static C_word
1963ilen(C_uword x) {
1964  C_uword y;
1965  C_word n = 0;
1966
1967#ifdef C_SIXTY_FOUR
1968  y = x >> 32; if (y != 0) { n += 32; x = y; }
1969#endif
1970  y = x >> 16; if (y != 0) { n += 16; x = y; }
1971  y = x >>  8; if (y != 0) { n +=  8; x = y; }
1972  y = x >>  4; if (y != 0) { n +=  4; x = y; }
1973  y = x >>  2; if (y != 0) { n +=  2; x = y; }
1974  y = x >>  1; if (y != 0) return n + 2;
1975  return n + x;
1976}
1977
1978static void
1979int_length(C_word c, C_word self, C_word k, C_word x)
1980{
1981  if (x & C_FIXNUM_BIT) {
1982    x = C_unfix(x);
1983    C_kontinue(k, C_fix(ilen((x < 0) ? ~x : x)));
1984  } else { /* bignum */
1985    bignum_type bigx;
1986    bignum_digit_type bd;
1987    int len_1;
1988
1989    bigx = big_of(x);
1990    if (BIGNUM_NEGATIVE_P(bigx)) {
1991      /* TODO: This can probably be done faster, by juggling with borrow */
1992      bignum_type notx = bignum_bitwise_not(bigx);
1993      len_1 = BIGNUM_LENGTH(notx) - 1;
1994      bd = *(BIGNUM_START_PTR(notx) + len_1); /* Most significant digit */
1995      BIGNUM_DEALLOCATE(notx);
1996    } else {
1997      len_1 = BIGNUM_LENGTH(bigx) - 1;
1998      bd = *(BIGNUM_START_PTR(bigx) + len_1); /* Most significant digit */
1999    }
2000    C_kontinue(k, C_fix(ilen(bd) + len_1 * BIGNUM_DIGIT_LENGTH));
2001  }
2002}
2003
2004static void 
2005int_shift_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
2006{
2007  bignum_type bigx, result;
2008
2009  if (y == C_fix(0)) C_kontinue(k, x); /* Done too (no shift) */
2010 
2011  /* Ensure x is a bignum */
2012  if((x & C_FIXNUM_BIT) != 0) {
2013    if (x == C_fix(0)) /* Skip everything else */
2014      C_kontinue(k, x);
2015   
2016    bigx = bignum_allocate_from_fixnum(x);
2017  } else {
2018    bigx = big_of(x);
2019  }
2020   
2021  result = bignum_arithmetic_shift(bigx, C_unfix(y));
2022  if ((x & C_FIXNUM_BIT) != 0)
2023    BIGNUM_DEALLOCATE(bigx);
2024  C_return_bignum(k, result);
2025}
2026
2027static bignum_type
2028bignum_arithmetic_shift(bignum_type arg1, C_word n)
2029{
2030  bignum_type tmp1, tmp2, result;
2031  if (BIGNUM_NEGATIVE_P(arg1) && n < 0) {
2032    tmp1 = bignum_bitwise_not(arg1);
2033    tmp2 = bignum_magnitude_ash(tmp1, n);
2034    BIGNUM_DEALLOCATE(tmp1);
2035    result = bignum_bitwise_not(tmp2);
2036    BIGNUM_DEALLOCATE(tmp2);
2037    return result;
2038  } else {
2039    return bignum_magnitude_ash(arg1, n);
2040  }
2041}
2042
2043static bignum_type
2044bignum_magnitude_ash(bignum_type arg1, C_word n)
2045{
2046  bignum_type result;
2047  bignum_digit_type *scan1;
2048  bignum_digit_type *scanr;
2049  bignum_digit_type *end;
2050
2051  C_word digit_offset,bit_offset;
2052  assert(n != 0); /* int_shift_fix checks this */
2053
2054  if (n > 0) {
2055    digit_offset = n / BIGNUM_DIGIT_LENGTH;
2056    bit_offset =   n % BIGNUM_DIGIT_LENGTH;
2057   
2058    result = bignum_allocate_zeroed (BIGNUM_LENGTH (arg1) + digit_offset + 1,
2059                                     BIGNUM_NEGATIVE_P(arg1));
2060
2061    scanr = BIGNUM_START_PTR (result) + digit_offset;
2062    scan1 = BIGNUM_START_PTR (arg1);
2063    end = scan1 + BIGNUM_LENGTH (arg1);
2064   
2065    while (scan1 < end) {
2066      *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset;
2067      *scanr = *scanr & BIGNUM_DIGIT_MASK;
2068      scanr++;
2069      *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset);
2070      *scanr = *scanr & BIGNUM_DIGIT_MASK;
2071    }
2072  }
2073  else if (n < 0
2074           && (-n >= (BIGNUM_LENGTH (arg1) * (bignum_length_type) BIGNUM_DIGIT_LENGTH)))
2075    result = BIGNUM_ZERO ();
2076
2077  else /* if (n < 0) */ {
2078    digit_offset = -n / BIGNUM_DIGIT_LENGTH;
2079    bit_offset =   -n % BIGNUM_DIGIT_LENGTH;
2080   
2081    result = bignum_allocate_zeroed (BIGNUM_LENGTH (arg1) - digit_offset,
2082                                     BIGNUM_NEGATIVE_P(arg1));
2083   
2084    scanr = BIGNUM_START_PTR (result);
2085    scan1 = BIGNUM_START_PTR (arg1) + digit_offset;
2086    end = scanr + BIGNUM_LENGTH (result) - 1;
2087   
2088    while (scanr < end) {
2089      *scanr =  (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset ;
2090      *scanr = (*scanr | 
2091        *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) & BIGNUM_DIGIT_MASK;
2092      scanr++;
2093    }
2094    *scanr =  (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset ;
2095  }
2096 
2097  return (bignum_trim (result));
2098}
2099
2100static bignum_type
2101bignum_pospos_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
2102{
2103  bignum_type result;
2104  bignum_length_type max_length;
2105
2106  bignum_digit_type *scan1, *end1, digit1;
2107  bignum_digit_type *scan2, *end2, digit2;
2108  bignum_digit_type *scanr, *endr;
2109
2110  max_length =  (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2))
2111               ? BIGNUM_LENGTH(arg1) : BIGNUM_LENGTH(arg2);
2112
2113  result = bignum_allocate(max_length, 0);
2114
2115  scanr = BIGNUM_START_PTR(result);
2116  scan1 = BIGNUM_START_PTR(arg1);
2117  scan2 = BIGNUM_START_PTR(arg2);
2118  endr = scanr + max_length;
2119  end1 = scan1 + BIGNUM_LENGTH(arg1);
2120  end2 = scan2 + BIGNUM_LENGTH(arg2);
2121
2122  while (scanr < endr) {
2123    digit1 = (scan1 < end1) ? *scan1++ : 0;
2124    digit2 = (scan2 < end2) ? *scan2++ : 0;
2125    /*
2126    fprintf(stderr, "[pospos op = %d, i = %ld, d1 = %lx, d2 = %lx]\n",
2127            op, endr - scanr, digit1, digit2);
2128            */
2129    *scanr++ = (op == bignum_and_op) ? digit1 & digit2 :
2130               (op == bignum_ior_op) ? digit1 | digit2 :
2131                                       digit1 ^ digit2;
2132  }
2133  return bignum_trim(result);
2134}
2135
2136static bignum_type
2137bignum_posneg_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
2138{
2139  bignum_type result;
2140  bignum_length_type max_length;
2141
2142  bignum_digit_type *scan1, *end1, digit1;
2143  bignum_digit_type *scan2, *end2, digit2, carry2;
2144  bignum_digit_type *scanr, *endr;
2145
2146  char neg_p = op == bignum_ior_op || op == bignum_xor_op;
2147
2148  max_length =  (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1)
2149               ? BIGNUM_LENGTH(arg1) : BIGNUM_LENGTH(arg2) + 1;
2150
2151  result = bignum_allocate(max_length, neg_p);
2152
2153  scanr = BIGNUM_START_PTR(result);
2154  scan1 = BIGNUM_START_PTR(arg1);
2155  scan2 = BIGNUM_START_PTR(arg2);
2156  endr = scanr + max_length;
2157  end1 = scan1 + BIGNUM_LENGTH(arg1);
2158  end2 = scan2 + BIGNUM_LENGTH(arg2);
2159
2160  carry2 = 1;
2161
2162  while (scanr < endr) {
2163    digit1 = (scan1 < end1) ? *scan1++ : 0;
2164    digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK)
2165             + carry2;
2166
2167    if (digit2 < BIGNUM_RADIX)
2168      carry2 = 0;
2169    else
2170      {
2171        digit2 = (digit2 - BIGNUM_RADIX);
2172        carry2 = 1;
2173      }
2174   
2175    *scanr++ = (op == bignum_and_op) ? digit1 & digit2 :
2176               (op == bignum_ior_op) ? digit1 | digit2 :
2177                                       digit1 ^ digit2;
2178  }
2179 
2180  if (neg_p)
2181    bignum_negate_magnitude(result);
2182
2183  return bignum_trim(result);
2184}
2185
2186static bignum_type
2187bignum_negneg_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
2188{
2189  bignum_type result;
2190  bignum_length_type max_length;
2191
2192  bignum_digit_type *scan1, *end1, digit1, carry1;
2193  bignum_digit_type *scan2, *end2, digit2, carry2;
2194  bignum_digit_type *scanr, *endr;
2195
2196  char neg_p = op == bignum_and_op || op == bignum_ior_op;
2197
2198  max_length =  (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2))
2199               ? BIGNUM_LENGTH(arg1) + 1 : BIGNUM_LENGTH(arg2) + 1;
2200
2201  result = bignum_allocate(max_length, neg_p);
2202
2203  scanr = BIGNUM_START_PTR(result);
2204  scan1 = BIGNUM_START_PTR(arg1);
2205  scan2 = BIGNUM_START_PTR(arg2);
2206  endr = scanr + max_length;
2207  end1 = scan1 + BIGNUM_LENGTH(arg1);
2208  end2 = scan2 + BIGNUM_LENGTH(arg2);
2209
2210  carry1 = 1;
2211  carry2 = 1;
2212
2213  while (scanr < endr) {
2214    digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1;
2215    digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
2216
2217    if (digit1 < BIGNUM_RADIX)
2218      carry1 = 0;
2219    else
2220      {
2221        digit1 = (digit1 - BIGNUM_RADIX);
2222        carry1 = 1;
2223      }
2224   
2225    if (digit2 < BIGNUM_RADIX)
2226      carry2 = 0;
2227    else
2228      {
2229        digit2 = (digit2 - BIGNUM_RADIX);
2230        carry2 = 1;
2231      }
2232   
2233    *scanr++ = (op == bignum_and_op) ? digit1 & digit2 :
2234               (op == bignum_ior_op) ? digit1 | digit2 :
2235                                       digit1 ^ digit2;
2236  }
2237
2238  if (neg_p)
2239    bignum_negate_magnitude(result);
2240
2241  return bignum_trim(result);
2242}
2243
2244static void
2245bignum_negate_magnitude(bignum_type arg)
2246{
2247  bignum_digit_type *scan;
2248  bignum_digit_type *end;
2249  bignum_digit_type digit;
2250  bignum_digit_type carry;
2251
2252  scan = BIGNUM_START_PTR(arg);
2253  end = scan + BIGNUM_LENGTH(arg);
2254
2255  carry = 1;
2256
2257  while (scan < end) {
2258    digit = (~*scan & BIGNUM_DIGIT_MASK) + carry;
2259
2260    if (digit < BIGNUM_RADIX)
2261      carry = 0;
2262    else
2263      {
2264        digit = (digit - BIGNUM_RADIX);
2265        carry = 1;
2266      }
2267   
2268    *scan++ = digit;
2269  }
2270}
2271
2272#define BIGNUM_STR_BLOCK_SIZE 8
2273
2274/*
2275 * Contains s48_bignum_to_digit_stream.  We may want to separate that out
2276 * again in case of huge numbers; it may be more efficient to write
2277 * those straight to files without going through the string conversion.
2278 */
2279static void
2280big_to_string(C_word c, C_word self, C_word k, C_word value, C_word radix)
2281{
2282  char *buf, *index;
2283  C_word len, counter=0;
2284  char *characters = "0123456789abcdef";
2285  C_word *tmp, ret;
2286  bignum_type bignum;
2287
2288  bignum = big_of(value);
2289  len = BIGNUM_STR_BLOCK_SIZE;
2290  buf = C_malloc(len);
2291  if (buf == NULL) {
2292    fprintf(stderr, "out of memory - can not allocate string");
2293    exit(EXIT_FAILURE);
2294  }
2295  index = buf + len - 1;
2296
2297  if (BIGNUM_NEGATIVE_P(bignum)) {
2298    *index-- = '-';
2299    counter++;
2300  }
2301
2302  radix = C_unfix(radix);
2303 
2304  assert((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
2305  if (! (BIGNUM_ZERO_P (bignum)))
2306    {
2307      bignum_type working_copy = (bignum_copy (bignum));
2308      bignum_digit_type * start = (BIGNUM_START_PTR (working_copy));
2309      bignum_digit_type * scan = (start + (BIGNUM_LENGTH (working_copy)));
2310      bignum_digit_type digit;
2311      while (start < scan)
2312        {
2313          if ((scan[-1]) == 0) {
2314            scan -= 1;
2315          } else {
2316            digit = bignum_destructive_scale_down (working_copy, radix);
2317            *index-- = characters[digit];
2318            if (++counter == len) {
2319              char *newbuf;
2320              newbuf = C_malloc(len + BIGNUM_STR_BLOCK_SIZE);
2321              if (newbuf == NULL) {
2322                fprintf(stderr, "out of memory - can not allocate string");
2323                exit(EXIT_FAILURE);
2324              }
2325              C_memcpy(newbuf + BIGNUM_STR_BLOCK_SIZE, buf, len);
2326              C_free(buf);
2327              buf = newbuf;
2328              index = newbuf + BIGNUM_STR_BLOCK_SIZE - 1;
2329              len += BIGNUM_STR_BLOCK_SIZE;
2330            }
2331          }
2332        }
2333      BIGNUM_DEALLOCATE (working_copy);
2334    }
2335
2336  if (BIGNUM_NEGATIVE_P(bignum))
2337    *index-- = '-';
2338
2339  tmp = C_alloc(C_SIZEOF_STRING(counter));
2340  ret = C_string(&tmp, counter, index+1);
2341  C_free(buf);
2342  C_kontinue(k, ret);
2343}
2344
2345static void
2346digits_to_big(C_word c, C_word self, C_word k, C_word n,
2347              C_word start, C_word end, C_word radix, C_word negp)
2348{
2349  char *str;
2350  size_t n_digits;
2351  int negative_p = (negp != C_SCHEME_FALSE);
2352  int digit;
2353  int hash = 0;
2354
2355  str = C_c_string(n) + C_unfix(start);
2356  n_digits = C_unfix(end)-C_unfix(start);
2357  radix = C_unfix(radix);
2358  hash = /* abs(radix / 2) */ 0;
2359 
2360#define DIGIT_TO_INT(x)         \
2361  (((x) == '#') ? hash :        \
2362   (((x) >= (int)'a') ?((x) - (int)'a' + 10) : ((x) - (int)'0')))
2363
2364  assert((radix > 1) && (radix < BIGNUM_RADIX_ROOT));
2365  if (n_digits == 0)
2366    C_kontinue(k, C_SCHEME_FALSE);
2367  if (n_digits == 1)
2368    {
2369      digit = DIGIT_TO_INT(C_tolower((int)*str));
2370      if (digit >= radix || digit < 0)
2371        C_kontinue(k, C_SCHEME_FALSE);
2372      else
2373        C_return_bignum(k, bignum_digit_to_bignum(digit, negative_p));
2374    }
2375  {
2376    bignum_length_type length;
2377    {
2378      unsigned int radix_copy = radix;
2379      unsigned int log_radix = 0;
2380      while (radix_copy > 0)
2381        {
2382          radix_copy >>= 1;
2383          log_radix += 1;
2384        }
2385      /* This length will be at least as large as needed. */
2386      length = (BIGNUM_BITS_TO_DIGITS (n_digits * log_radix));
2387    }
2388    {
2389      bignum_type result = (bignum_allocate_zeroed (length, negative_p));
2390      while ((n_digits--) > 0)
2391        {
2392          digit = DIGIT_TO_INT(C_tolower((int)*str));
2393          str++;
2394          if (digit >= radix || digit < 0) {
2395            BIGNUM_DEALLOCATE(result);
2396            C_kontinue(k, C_SCHEME_FALSE);
2397          }
2398       
2399          bignum_destructive_scale_up (result, ((bignum_digit_type) radix));
2400          bignum_destructive_add (result, digit);
2401        }
2402      C_return_bignum (k, bignum_trim (result));
2403    }
2404  }
2405}
Note: See TracBrowser for help on using the repository browser.