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

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

numbers: simplify and inline shorten_bignum into bignum_add_unsigned, then remove it; it was only used in one place and a bit silly because we always know the target bignum is one smaller, because the allocation allowed for a carry

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