1 /* -*-C-*-
2
3 Copyright (C) 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994,
4 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005,
5 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014 Massachusetts
6 Institute of Technology
7
8 This file is part of MIT/GNU Scheme.
9
10 MIT/GNU Scheme is free software; you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation; either version 2 of the License, or (at
13 your option) any later version.
14
15 MIT/GNU Scheme is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with MIT/GNU Scheme; if not, write to the Free Software
22 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301,
23 USA.
24
25 */
26
27 /* Implementation of Bignums (unlimited precision integers) */
28
29 #ifdef MIT_SCHEME
30 # include "scheme.h"
31 #else
32 # include "bignum.h"
33 #endif
34
35 #include "bignmint.h"
36 #include "bits.h"
37
38 #ifndef MIT_SCHEME
39
40 static bignum_type
bignum_malloc(bignum_length_type length)41 bignum_malloc (bignum_length_type length)
42 {
43 extern char * malloc ();
44 char * result = (malloc ((length + 1) * (sizeof (bignum_digit_type))));
45 BIGNUM_ASSERT (result != ((char *) 0));
46 return ((bignum_type) result);
47 }
48
49 static bignum_type
bignum_realloc(bignum_type bignum,bignum_length_type length)50 bignum_realloc (bignum_type bignum, bignum_length_type length)
51 {
52 extern char * realloc ();
53 char * result =
54 (realloc (((char *) bignum),
55 ((length + 1) * (sizeof (bignum_digit_type)))));
56 BIGNUM_ASSERT (result != ((char *) 0));
57 return ((bignum_type) result);
58 }
59
60 #endif /* not MIT_SCHEME */
61
62 /* Forward references */
63 static int bignum_equal_p_unsigned (bignum_type, bignum_type);
64 static enum bignum_comparison bignum_compare_unsigned
65 (bignum_type, bignum_type);
66 static bignum_type bignum_add_unsigned (bignum_type, bignum_type, int);
67 static bignum_type bignum_subtract_unsigned (bignum_type, bignum_type);
68 static bignum_type bignum_multiply_unsigned (bignum_type, bignum_type, int);
69 static bignum_type bignum_multiply_unsigned_small_factor
70 (bignum_type, bignum_digit_type, int);
71 static void bignum_destructive_scale_up (bignum_type, bignum_digit_type);
72 static void bignum_destructive_add (bignum_type, bignum_digit_type);
73 static void bignum_divide_unsigned_large_denominator
74 (bignum_type, bignum_type, bignum_type *, bignum_type *, int, int);
75 static void bignum_destructive_normalization (bignum_type, bignum_type, int);
76 static void bignum_destructive_unnormalization (bignum_type, int);
77 static void bignum_divide_unsigned_normalized
78 (bignum_type, bignum_type, bignum_type);
79 static bignum_digit_type bignum_divide_subtract
80 (bignum_digit_type *, bignum_digit_type *,
81 bignum_digit_type, bignum_digit_type *);
82 static void bignum_divide_unsigned_medium_denominator
83 (bignum_type, bignum_digit_type, bignum_type *, bignum_type *, int, int);
84 static bignum_digit_type bignum_digit_divide
85 (bignum_digit_type, bignum_digit_type,
86 bignum_digit_type, bignum_digit_type *);
87 static bignum_digit_type bignum_digit_divide_subtract
88 (bignum_digit_type, bignum_digit_type,
89 bignum_digit_type, bignum_digit_type *);
90 static void bignum_divide_unsigned_small_denominator
91 (bignum_type, bignum_digit_type, bignum_type *, bignum_type *, int, int);
92 static bignum_digit_type bignum_destructive_scale_down
93 (bignum_type, bignum_digit_type);
94 static bignum_type bignum_remainder_unsigned_small_denominator
95 (bignum_type, bignum_digit_type, int);
96 static bignum_type bignum_digit_to_bignum (bignum_digit_type, int);
97 static bignum_type bignum_allocate (bignum_length_type, int);
98 static bignum_type bignum_allocate_zeroed (bignum_length_type, int);
99 static bignum_type bignum_shorten_length (bignum_type, bignum_length_type);
100 static bignum_type bignum_trim (bignum_type);
101 static bignum_type bignum_copy (bignum_type);
102 static bignum_type bignum_new_sign (bignum_type, int);
103 static bignum_type bignum_maybe_new_sign (bignum_type, int);
104 static void bignum_destructive_copy (bignum_type, bignum_type);
105
106 /* Exports */
107
108 bignum_type
bignum_make_zero(void)109 bignum_make_zero (void)
110 {
111 bignum_type result = (BIGNUM_ALLOCATE (0));
112 BIGNUM_SET_HEADER (result, 0, 0);
113 return (result);
114 }
115
116 bignum_type
bignum_make_one(int negative_p)117 bignum_make_one (int negative_p)
118 {
119 bignum_type result = (BIGNUM_ALLOCATE (1));
120 BIGNUM_SET_HEADER (result, 1, negative_p);
121 (BIGNUM_REF (result, 0)) = 1;
122 return (result);
123 }
124
125 int
bignum_equal_p(bignum_type x,bignum_type y)126 bignum_equal_p (bignum_type x, bignum_type y)
127 {
128 return
129 ((BIGNUM_ZERO_P (x))
130 ? (BIGNUM_ZERO_P (y))
131 : ((! (BIGNUM_ZERO_P (y)))
132 && ((BIGNUM_NEGATIVE_P (x))
133 ? (BIGNUM_NEGATIVE_P (y))
134 : (! (BIGNUM_NEGATIVE_P (y))))
135 && (bignum_equal_p_unsigned (x, y))));
136 }
137
138 enum bignum_comparison
bignum_test(bignum_type bignum)139 bignum_test (bignum_type bignum)
140 {
141 return
142 ((BIGNUM_ZERO_P (bignum))
143 ? bignum_comparison_equal
144 : (BIGNUM_NEGATIVE_P (bignum))
145 ? bignum_comparison_less
146 : bignum_comparison_greater);
147 }
148
149 enum bignum_comparison
bignum_compare(bignum_type x,bignum_type y)150 bignum_compare (bignum_type x, bignum_type y)
151 {
152 return
153 ((BIGNUM_ZERO_P (x))
154 ? ((BIGNUM_ZERO_P (y))
155 ? bignum_comparison_equal
156 : (BIGNUM_NEGATIVE_P (y))
157 ? bignum_comparison_greater
158 : bignum_comparison_less)
159 : (BIGNUM_ZERO_P (y))
160 ? ((BIGNUM_NEGATIVE_P (x))
161 ? bignum_comparison_less
162 : bignum_comparison_greater)
163 : (BIGNUM_NEGATIVE_P (x))
164 ? ((BIGNUM_NEGATIVE_P (y))
165 ? (bignum_compare_unsigned (y, x))
166 : (bignum_comparison_less))
167 : ((BIGNUM_NEGATIVE_P (y))
168 ? (bignum_comparison_greater)
169 : (bignum_compare_unsigned (x, y))));
170 }
171
172 bignum_type
bignum_add(bignum_type x,bignum_type y)173 bignum_add (bignum_type x, bignum_type y)
174 {
175 return
176 ((BIGNUM_ZERO_P (x))
177 ? (BIGNUM_MAYBE_COPY (y))
178 : (BIGNUM_ZERO_P (y))
179 ? (BIGNUM_MAYBE_COPY (x))
180 : ((BIGNUM_NEGATIVE_P (x))
181 ? ((BIGNUM_NEGATIVE_P (y))
182 ? (bignum_add_unsigned (x, y, 1))
183 : (bignum_subtract_unsigned (y, x)))
184 : ((BIGNUM_NEGATIVE_P (y))
185 ? (bignum_subtract_unsigned (x, y))
186 : (bignum_add_unsigned (x, y, 0)))));
187 }
188
189 bignum_type
bignum_subtract(bignum_type x,bignum_type y)190 bignum_subtract (bignum_type x, bignum_type y)
191 {
192 return
193 ((BIGNUM_ZERO_P (x))
194 ? ((BIGNUM_ZERO_P (y))
195 ? (BIGNUM_MAYBE_COPY (y))
196 : (bignum_new_sign (y, (! (BIGNUM_NEGATIVE_P (y))))))
197 : ((BIGNUM_ZERO_P (y))
198 ? (BIGNUM_MAYBE_COPY (x))
199 : ((BIGNUM_NEGATIVE_P (x))
200 ? ((BIGNUM_NEGATIVE_P (y))
201 ? (bignum_subtract_unsigned (y, x))
202 : (bignum_add_unsigned (x, y, 1)))
203 : ((BIGNUM_NEGATIVE_P (y))
204 ? (bignum_add_unsigned (x, y, 0))
205 : (bignum_subtract_unsigned (x, y))))));
206 }
207
208 bignum_type
bignum_negate(bignum_type x)209 bignum_negate (bignum_type x)
210 {
211 return
212 ((BIGNUM_ZERO_P (x))
213 ? (BIGNUM_MAYBE_COPY (x))
214 : (bignum_new_sign (x, (! (BIGNUM_NEGATIVE_P (x))))));
215 }
216
217 bignum_type
bignum_multiply(bignum_type x,bignum_type y)218 bignum_multiply (bignum_type x, bignum_type y)
219 {
220 bignum_length_type x_length = (BIGNUM_LENGTH (x));
221 bignum_length_type y_length = (BIGNUM_LENGTH (y));
222 int negative_p =
223 ((BIGNUM_NEGATIVE_P (x))
224 ? (! (BIGNUM_NEGATIVE_P (y)))
225 : (BIGNUM_NEGATIVE_P (y)));
226 if (BIGNUM_ZERO_P (x))
227 return (BIGNUM_MAYBE_COPY (x));
228 if (BIGNUM_ZERO_P (y))
229 return (BIGNUM_MAYBE_COPY (y));
230 if (x_length == 1)
231 {
232 bignum_digit_type digit = (BIGNUM_REF (x, 0));
233 if (digit == 1)
234 return (bignum_maybe_new_sign (y, negative_p));
235 if (digit < BIGNUM_RADIX_ROOT)
236 return (bignum_multiply_unsigned_small_factor (y, digit, negative_p));
237 }
238 if (y_length == 1)
239 {
240 bignum_digit_type digit = (BIGNUM_REF (y, 0));
241 if (digit == 1)
242 return (bignum_maybe_new_sign (x, negative_p));
243 if (digit < BIGNUM_RADIX_ROOT)
244 return (bignum_multiply_unsigned_small_factor (x, digit, negative_p));
245 }
246 return (bignum_multiply_unsigned (x, y, negative_p));
247 }
248
249 int
bignum_divide(bignum_type numerator,bignum_type denominator,bignum_type * quotient,bignum_type * remainder)250 bignum_divide (bignum_type numerator, bignum_type denominator,
251 bignum_type * quotient, bignum_type * remainder)
252 {
253 if (BIGNUM_ZERO_P (denominator))
254 return (1);
255 if (BIGNUM_ZERO_P (numerator))
256 {
257 (*quotient) = (BIGNUM_MAYBE_COPY (numerator));
258 (*remainder) = (BIGNUM_MAYBE_COPY (numerator));
259 }
260 else
261 {
262 int r_negative_p = (BIGNUM_NEGATIVE_P (numerator));
263 int q_negative_p =
264 ((BIGNUM_NEGATIVE_P (denominator)) ? (! r_negative_p) : r_negative_p);
265 switch (bignum_compare_unsigned (numerator, denominator))
266 {
267 case bignum_comparison_equal:
268 {
269 (*quotient) = (BIGNUM_ONE (q_negative_p));
270 (*remainder) = (BIGNUM_ZERO ());
271 break;
272 }
273 case bignum_comparison_less:
274 {
275 (*quotient) = (BIGNUM_ZERO ());
276 (*remainder) = (BIGNUM_MAYBE_COPY (numerator));
277 break;
278 }
279 case bignum_comparison_greater:
280 {
281 if ((BIGNUM_LENGTH (denominator)) == 1)
282 {
283 bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
284 if (digit == 1)
285 {
286 (*quotient) =
287 (bignum_maybe_new_sign (numerator, q_negative_p));
288 (*remainder) = (BIGNUM_ZERO ());
289 break;
290 }
291 else if (digit < BIGNUM_RADIX_ROOT)
292 {
293 bignum_divide_unsigned_small_denominator
294 (numerator, digit,
295 quotient, remainder,
296 q_negative_p, r_negative_p);
297 break;
298 }
299 else
300 {
301 bignum_divide_unsigned_medium_denominator
302 (numerator, digit,
303 quotient, remainder,
304 q_negative_p, r_negative_p);
305 break;
306 }
307 }
308 bignum_divide_unsigned_large_denominator
309 (numerator, denominator,
310 quotient, remainder,
311 q_negative_p, r_negative_p);
312 break;
313 }
314 }
315 }
316 return (0);
317 }
318
319 bignum_type
bignum_quotient(bignum_type numerator,bignum_type denominator)320 bignum_quotient (bignum_type numerator, bignum_type denominator)
321 {
322 if (BIGNUM_ZERO_P (denominator))
323 return (BIGNUM_OUT_OF_BAND);
324 if (BIGNUM_ZERO_P (numerator))
325 return (BIGNUM_MAYBE_COPY (numerator));
326 {
327 int q_negative_p =
328 ((BIGNUM_NEGATIVE_P (denominator))
329 ? (! (BIGNUM_NEGATIVE_P (numerator)))
330 : (BIGNUM_NEGATIVE_P (numerator)));
331 switch (bignum_compare_unsigned (numerator, denominator))
332 {
333 case bignum_comparison_equal:
334 return (BIGNUM_ONE (q_negative_p));
335 case bignum_comparison_less:
336 return (BIGNUM_ZERO ());
337 case bignum_comparison_greater:
338 {
339 bignum_type quotient;
340 if ((BIGNUM_LENGTH (denominator)) == 1)
341 {
342 bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
343 if (digit == 1)
344 return (bignum_maybe_new_sign (numerator, q_negative_p));
345 if (digit < BIGNUM_RADIX_ROOT)
346 bignum_divide_unsigned_small_denominator
347 (numerator, digit,
348 ("ient), ((bignum_type *) 0),
349 q_negative_p, 0);
350 else
351 bignum_divide_unsigned_medium_denominator
352 (numerator, digit,
353 ("ient), ((bignum_type *) 0),
354 q_negative_p, 0);
355 }
356 else
357 bignum_divide_unsigned_large_denominator
358 (numerator, denominator,
359 ("ient), ((bignum_type *) 0),
360 q_negative_p, 0);
361 return (quotient);
362 }
363 default:
364 /*NOTREACHED*/
365 return (0);
366 }
367 }
368 }
369
370 bignum_type
bignum_remainder(bignum_type numerator,bignum_type denominator)371 bignum_remainder (bignum_type numerator, bignum_type denominator)
372 {
373 if (BIGNUM_ZERO_P (denominator))
374 return (BIGNUM_OUT_OF_BAND);
375 if (BIGNUM_ZERO_P (numerator))
376 return (BIGNUM_MAYBE_COPY (numerator));
377 switch (bignum_compare_unsigned (numerator, denominator))
378 {
379 case bignum_comparison_equal:
380 return (BIGNUM_ZERO ());
381 case bignum_comparison_less:
382 return (BIGNUM_MAYBE_COPY (numerator));
383 case bignum_comparison_greater:
384 {
385 bignum_type remainder;
386 if ((BIGNUM_LENGTH (denominator)) == 1)
387 {
388 bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
389 if ((digit & (digit-1)) == 0) /* i.e. digit = 2^k, including 1 */
390 {
391 bignum_digit_type unsigned_remainder;
392 if (BIGNUM_LENGTH (numerator) == 0)
393 return (BIGNUM_ZERO ());
394 unsigned_remainder = (digit-1) & (BIGNUM_REF (numerator, 0));
395 if (unsigned_remainder == 0)
396 return (BIGNUM_ZERO ());
397 return (bignum_digit_to_bignum
398 (unsigned_remainder, BIGNUM_NEGATIVE_P (numerator)));
399 }
400 if (digit < BIGNUM_RADIX_ROOT)
401 return
402 (bignum_remainder_unsigned_small_denominator
403 (numerator, digit, (BIGNUM_NEGATIVE_P (numerator))));
404 bignum_divide_unsigned_medium_denominator
405 (numerator, digit,
406 ((bignum_type *) 0), (&remainder),
407 0, (BIGNUM_NEGATIVE_P (numerator)));
408 }
409 else
410 bignum_divide_unsigned_large_denominator
411 (numerator, denominator,
412 ((bignum_type *) 0), (&remainder),
413 0, (BIGNUM_NEGATIVE_P (numerator)));
414 return (remainder);
415 }
416 default:
417 /*NOTREACHED*/
418 return (0);
419 }
420 }
421
422 static bignum_type
bignum_from_digits(bignum_digit_type * start_digits,bignum_digit_type * end_digits,bool negative_p)423 bignum_from_digits (bignum_digit_type *start_digits,
424 bignum_digit_type *end_digits,
425 bool negative_p)
426 {
427 bignum_type result =
428 (bignum_allocate ((end_digits - start_digits), negative_p));
429 bignum_digit_type *scan_digits = start_digits;
430 bignum_digit_type *scan_result = (BIGNUM_START_PTR (result));
431 while (scan_digits < end_digits)
432 (*scan_result++) = (*scan_digits++);
433 return (result);
434 }
435
436 #define DEFINE_INT_TO_BIGNUM(NAME, TYPE, UTYPE) \
437 bignum_type \
438 NAME (TYPE n) \
439 { \
440 /* Special cases win when these small constants are cached. */ \
441 if (n == 0) return (BIGNUM_ZERO ()); \
442 if (n == 1) return (BIGNUM_ONE (0)); \
443 if (n == -1) return (BIGNUM_ONE (1)); \
444 { \
445 int negative_p; \
446 bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_TYPE (TYPE)]; \
447 bignum_digit_type * end_digits = result_digits; \
448 UTYPE accumulator = ((negative_p = (n < 0)) ? (-n) : n); \
449 do { \
450 (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK); \
451 accumulator >>= BIGNUM_DIGIT_LENGTH; \
452 } while (accumulator != 0); \
453 return \
454 (bignum_from_digits (result_digits, end_digits, negative_p)); \
455 } \
456 }
457
DEFINE_INT_TO_BIGNUM(long_to_bignum,long,unsigned long)458 DEFINE_INT_TO_BIGNUM (long_to_bignum, long, unsigned long)
459 DEFINE_INT_TO_BIGNUM (intmax_to_bignum, intmax_t, uintmax_t)
460
461 #define DEFINE_BIGNUM_TO_INT(NAME, TYPE, UTYPE) \
462 TYPE \
463 NAME (bignum_type bignum) \
464 { \
465 if (BIGNUM_ZERO_P (bignum)) \
466 return (0); \
467 { \
468 UTYPE accumulator = 0; \
469 bignum_digit_type * start = (BIGNUM_START_PTR (bignum)); \
470 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum))); \
471 while (start < scan) \
472 accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan)); \
473 return \
474 ((BIGNUM_NEGATIVE_P (bignum)) \
475 ? (- ((TYPE) accumulator)) \
476 : ((TYPE) accumulator)); \
477 } \
478 }
479
480 DEFINE_BIGNUM_TO_INT (bignum_to_long, long, unsigned long)
481 DEFINE_BIGNUM_TO_INT (bignum_to_intmax, intmax_t, uintmax_t)
482
483 #define DEFINE_UINT_TO_BIGNUM(NAME, TYPE) \
484 bignum_type \
485 NAME (TYPE n) \
486 { \
487 /* Special cases win when these small constants are cached. */ \
488 if (n == 0) return (BIGNUM_ZERO ()); \
489 if (n == 1) return (BIGNUM_ONE (0)); \
490 { \
491 bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_TYPE (TYPE)]; \
492 bignum_digit_type * end_digits = result_digits; \
493 TYPE accumulator = n; \
494 do { \
495 (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK); \
496 accumulator >>= BIGNUM_DIGIT_LENGTH; \
497 } while (accumulator != 0); \
498 return (bignum_from_digits (result_digits, end_digits, false)); \
499 } \
500 }
501
502 DEFINE_UINT_TO_BIGNUM (ulong_to_bignum, unsigned long)
503 DEFINE_UINT_TO_BIGNUM (uintmax_to_bignum, uintmax_t)
504
505 #define DEFINE_BIGNUM_TO_UINT(NAME, TYPE) \
506 TYPE \
507 NAME (bignum_type bignum) \
508 { \
509 if (BIGNUM_ZERO_P (bignum)) \
510 return (0); \
511 BIGNUM_ASSERT (BIGNUM_POSITIVE_P (bignum)); \
512 { \
513 TYPE accumulator = 0; \
514 bignum_digit_type * start = (BIGNUM_START_PTR (bignum)); \
515 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum))); \
516 while (start < scan) \
517 accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan)); \
518 return (accumulator); \
519 } \
520 }
521
522 DEFINE_BIGNUM_TO_UINT (bignum_to_ulong, unsigned long)
523 DEFINE_BIGNUM_TO_UINT (bignum_to_uintmax, uintmax_t)
524
525 #define DTB_WRITE_DIGIT(n_bits) do \
526 { \
527 significand *= (1UL << (n_bits)); \
528 digit = ((bignum_digit_type) significand); \
529 (*--scan) = digit; \
530 significand -= ((double) digit); \
531 n_valid_bits -= (n_bits); \
532 } while (0)
533
534 bignum_type
535 double_to_bignum (double x)
536 {
537 int exponent;
538 double significand = (frexp (x, (&exponent)));
539 if (exponent <= 0) return (BIGNUM_ZERO ());
540 if (exponent == 1) return (BIGNUM_ONE (x < 0));
541 if (significand < 0) significand = (-significand);
542 {
543 bignum_length_type length = (BIGNUM_BITS_TO_DIGITS (exponent));
544 bignum_type result = (bignum_allocate (length, (x < 0)));
545 bignum_digit_type * start = (BIGNUM_START_PTR (result));
546 bignum_digit_type * scan = (start + length);
547 unsigned int n_valid_bits = DBL_MANT_DIG;
548 bignum_digit_type digit;
549 {
550 int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
551 if (odd_bits > 0)
552 DTB_WRITE_DIGIT (odd_bits);
553 }
554 while (start < scan)
555 {
556 if ((significand == 0) || (n_valid_bits == 0))
557 {
558 while (start < scan)
559 (*--scan) = 0;
560 break;
561 }
562 if (n_valid_bits >= BIGNUM_DIGIT_LENGTH)
563 DTB_WRITE_DIGIT (BIGNUM_DIGIT_LENGTH);
564 else
565 {
566 significand *= (1UL << BIGNUM_DIGIT_LENGTH);
567 digit = ((bignum_digit_type) significand);
568 (*--scan)
569 = (digit
570 & (((1UL << n_valid_bits) - 1UL)
571 << (BIGNUM_DIGIT_LENGTH - n_valid_bits)));
572 significand = 0.0;
573 n_valid_bits = 0;
574 }
575 }
576 return (result);
577 }
578 }
579
580 #undef DTB_WRITE_DIGIT
581
582 /* This version sometimes gets the answer wrong due to rounding errors.
583 It would be slightly better if the digits were accumulated lsb to msb.
584 */
585 /*
586 double
587 bignum_to_double (bignum_type bignum)
588 {
589 if (BIGNUM_ZERO_P (bignum))
590 return (0);
591 {
592 double accumulator = 0;
593 bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
594 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
595 while (start < scan)
596 accumulator = ((accumulator * BIGNUM_RADIX) + (*--scan));
597 return ((BIGNUM_NEGATIVE_P (bignum)) ? (-accumulator) : accumulator);
598 }
599 }
600 */
601
602 double
bignum_to_double(bignum_type bignum)603 bignum_to_double (bignum_type bignum)
604 {
605 if (BIGNUM_ZERO_P (bignum))
606 return (0.0);
607
608 {
609 bignum_length_type length = (BIGNUM_LENGTH (bignum));
610 bignum_length_type index = length - 1;
611 bignum_length_type scale_words;
612 bignum_digit_type msd = (BIGNUM_REF (bignum, (index)));
613 #if (FLT_RADIX == 2)
614 int bits_to_get = DBL_MANT_DIG; /* includes implicit 1 */
615 #else
616 # include "error: must have FLT_RADIX==2"
617 #endif
618 double value = 0;
619 bignum_digit_type guard_bit_mask = BIGNUM_RADIX>>1;
620 bignum_digit_type rounding_correction = 0;
621 int current_digit_bit_count = (ulong_length_in_bits (msd));
622 bignum_digit_type mask = (BIGNUM_DIGIT_ONES (current_digit_bit_count));
623
624 while (1) {
625 if (current_digit_bit_count > bits_to_get) {
626 guard_bit_mask = (1UL << (current_digit_bit_count - bits_to_get - 1));
627 mask &= ~((guard_bit_mask << 1) - 1UL);
628 current_digit_bit_count = bits_to_get;
629 bits_to_get = 0;
630 } else {
631 bits_to_get -= current_digit_bit_count;
632 }
633
634 value = (value * BIGNUM_RADIX) + ((BIGNUM_REF (bignum, index)) & mask);
635
636 if (bits_to_get == 0) {
637 scale_words = index;
638 if (current_digit_bit_count == BIGNUM_DIGIT_LENGTH) {
639 if (index == 0) /* there is no guard bit */
640 goto finished;
641 guard_bit_mask = (1UL << (BIGNUM_DIGIT_LENGTH - 1));
642 rounding_correction = 1;
643 index -= 1;
644 } else {
645 rounding_correction = (guard_bit_mask << 1);
646 }
647 break;
648 }
649 if (index == 0) /* fewer than DBL_MANT_DIG bits */
650 goto finished;
651
652 index -= 1;
653 current_digit_bit_count = BIGNUM_DIGIT_LENGTH;
654 mask = BIGNUM_DIGIT_MASK;
655 }
656
657 /* round-to-even depending on lsb, guard and following bits: lgfffff */
658
659 if ((BIGNUM_REF(bignum,index) & guard_bit_mask) == 0) /* case x0xxxx */
660 goto round_down;
661
662 if ((BIGNUM_REF(bignum,index) & (guard_bit_mask-1)) != 0) /* case x1xx1x */
663 goto round_up;
664
665 /* cases 110000 and 1101xx: test "odd?", i.e. round-to-even rounds up */
666 if ((guard_bit_mask << 1) == BIGNUM_RADIX) {
667 if (((BIGNUM_REF (bignum, index+1)) & 1UL) != 0) /* "odd?" */
668 goto round_up;
669 } else {
670 if (((BIGNUM_REF (bignum, index)) & (guard_bit_mask << 1)) != 0)
671 goto round_up;
672 }
673
674 if (index==0) /* case 010000, no more words of following bits */
675 goto finished;
676
677 { /* distinguish between cases 0100...00 and 0100..1xx, multiple words */
678 bignum_length_type index2 = index - 1;
679 while (index2 >= 0) {
680 if ((BIGNUM_REF (bignum, index2)) != 0)
681 goto round_up;
682 index2--;
683 }
684 goto round_down;
685 }
686
687 round_up:
688 value += rounding_correction;
689 round_down:
690 /* note, ldexp `sticks' at the maximal non-infinity value, which
691 is a reasonable compromise for numbers with DBL_MAX_EXP bits
692 that round up */
693 if (scale_words > 0)
694 value = ldexp (value, scale_words * BIGNUM_DIGIT_LENGTH);
695
696 finished:
697 return ((BIGNUM_NEGATIVE_P (bignum)) ? (-value) : value);
698 }
699 }
700
701 /*
702 ;; Code to test bignum_to_double
703
704 (declare (usual-integrations))
705
706 (define integer->flonum (make-primitive-procedure 'integer->flonum 2))
707
708 (define (check n)
709 (let ((f1 (integer->flonum n #b10))
710 (f2 (exact->inexact n)))
711 (if (and f1 (= f1 f2))
712 (number->string f1 2)
713 (begin
714 (pp n)
715 (pp (number->string f1 2))
716 (pp (number->string f2 2))))))
717
718 (define (test)
719 (define n 0)
720 (do ((i 0 (+ i 1))) ; guard bit zone patterns
721 ((= i 256))
722 (do ((j 0 (+ j 1))) ; general word alignment
723 ((= j 30))
724 (do ((e 0 (+ e 1))) ; random `insignificant' bit
725 ((= e 2))
726 (do ((up 0 (+ up 1))) ; test potential for carry
727 ((= up 2))
728 (do ((end-pad 0 (+ end-pad 100)))
729 ((> end-pad 100))
730 (set! n (+ n 1)) (if (= (remainder n 1000) 0) (pp `(test ,n)))
731 (let ((s (string-append "1"
732 (make-string 48 (vector-ref '#(#\0 #\1) up))
733 (string-pad-left (number->string i 2) 8 #\0)
734 (make-string (* j 23) #\0) ; gcd(23,30)=1
735 (number->string e)
736 (make-string end-pad #\0))))
737 (check (string->number s 2)))))))))
738 */
739
740 int
bignum_fits_in_word_p(bignum_type bignum,long word_length,int twos_complement_p)741 bignum_fits_in_word_p (bignum_type bignum, long word_length,
742 int twos_complement_p)
743 {
744 unsigned int n_bits;
745
746 if (twos_complement_p)
747 n_bits = (word_length - 1);
748 else
749 {
750 if (BIGNUM_NEGATIVE_P (bignum))
751 return (0);
752 n_bits = word_length;
753 }
754
755 BIGNUM_ASSERT (n_bits > 0);
756 {
757 bignum_length_type length = (BIGNUM_LENGTH (bignum));
758 unsigned int max_digits = (BIGNUM_BITS_TO_DIGITS (n_bits));
759 if (((unsigned int) length) < max_digits)
760 return (1);
761 if (((unsigned int) length) > max_digits)
762 return (0);
763 {
764 bignum_digit_type msd = (BIGNUM_REF (bignum, (length - 1)));
765 bignum_digit_type max
766 = (1UL << (n_bits - ((length - 1) * BIGNUM_DIGIT_LENGTH)));
767 return
768 (((msd < max)
769 || (twos_complement_p
770 && (msd == max)
771 && (BIGNUM_NEGATIVE_P (bignum)))));
772 }
773 }
774 }
775
776 /* All these bitwise operations are complicated because they interpret
777 integers as their two's-complement representations, whereas bignums
778 are stored in a ones-complement representation. */
779
780 bignum_type
bignum_bitwise_not(bignum_type bignum)781 bignum_bitwise_not (bignum_type bignum)
782 {
783 return (bignum_subtract ((BIGNUM_ONE (1)), bignum));
784 }
785
786 #define DEFINE_BITWISE_UNSIGNED(NAME, COMMUTE, OP, POSTOP) \
787 static bignum_type \
788 NAME (bignum_type x, bignum_type y) \
789 { \
790 if ((BIGNUM_LENGTH (x)) < (BIGNUM_LENGTH (y))) \
791 COMMUTE (y, x); \
792 { \
793 bignum_length_type x_length = (BIGNUM_LENGTH (x)); \
794 bignum_length_type y_length = (BIGNUM_LENGTH (y)); \
795 bignum_length_type r_length = x_length; \
796 bignum_type r = (bignum_allocate (r_length, 0)); \
797 bignum_digit_type *x_scan = (BIGNUM_START_PTR (x)); \
798 bignum_digit_type *x_end = (x_scan + x_length); \
799 bignum_digit_type *y_scan = (BIGNUM_START_PTR (y)); \
800 bignum_digit_type *y_end = (y_scan + y_length); \
801 bignum_digit_type *r_scan = (BIGNUM_START_PTR (r)); \
802 BIGNUM_ASSERT (x_length >= y_length); \
803 while (y_scan < y_end) \
804 (*r_scan++) = (BITWISE_##OP ((*x_scan++), (*y_scan++))); \
805 while (x_scan < x_end) \
806 (*r_scan++) = (BITWISE_##OP ((*x_scan++), 0)); \
807 return (POSTOP (bignum_trim (r))); \
808 } \
809 }
810
811 #define BITWISE_AND(x, y) ((x) & (y))
812 #define BITWISE_ANDC2(x, y) ((x) &~ (y))
813 #define BITWISE_ANDC1(x, y) ((y) &~ (x))
814 #define BITWISE_XOR(x, y) ((x) ^ (y))
815 #define BITWISE_IOR(x, y) ((x) | (y))
816
817 /* These don't work, because they set the high bits. */
818 /* #define BITWISE_ORC2(x, y) (BIGNUM_DIGIT_MASK & ((x) |~ (y))) */
819 /* #define BITWISE_ORC1(x, y) (BIGNUM_DIGIT_MASK & ((y) |~ (x))) */
820
821 /* Kludgey syntactic hack! */
822 #define COMMUTE(name) return bignum_bitwise_##name##_unsigned
823 #define SWAP(x, y) do { bignum_type t = x; x = y; y = t; } while (0)
824
825 static bignum_type bignum_bitwise_andc1_unsigned (bignum_type, bignum_type);
826 static bignum_type bignum_bitwise_orc1_unsigned (bignum_type, bignum_type);
827
828 /* These definitions are ordered by their truth tables. */
829
830 /* DEFINE_BITWISE_UNSIGNED (bignum_bitwise_clear_unsigned, SWAP, CLEAR,) */
831 DEFINE_BITWISE_UNSIGNED (bignum_bitwise_and_unsigned, SWAP, AND,)
832 DEFINE_BITWISE_UNSIGNED (bignum_bitwise_andc2_unsigned, COMMUTE(andc1), ANDC2,)
833 /* DEFINE_BITWISE_UNSIGNED (bignum_bitwise_arg1_unsigned, COMMUTE(ARG2), ARG1,) */
834 DEFINE_BITWISE_UNSIGNED (bignum_bitwise_andc1_unsigned, COMMUTE(andc2), ANDC1,)
835 /* DEFINE_BITWISE_UNSIGNED (bignum_bitwise_arg2_unsigned, ARG2,) */
836 DEFINE_BITWISE_UNSIGNED (bignum_bitwise_xor_unsigned, SWAP, XOR,)
837 DEFINE_BITWISE_UNSIGNED (bignum_bitwise_ior_unsigned, SWAP, IOR,)
DEFINE_BITWISE_UNSIGNED(bignum_bitwise_nor_unsigned,SWAP,IOR,bignum_bitwise_not)838 DEFINE_BITWISE_UNSIGNED (bignum_bitwise_nor_unsigned, SWAP, IOR, bignum_bitwise_not)
839 DEFINE_BITWISE_UNSIGNED (bignum_bitwise_eqv_unsigned, SWAP, XOR, bignum_bitwise_not)
840 /* DEFINE_BITWISE_UNSIGNED (bignum_bitwise_not2_unsigned, COMMUTE(not1), ARG1, bignum_bitwise_not) */
841 /* DEFINE_BITWISE_UNSIGNED (bignum_bitwise_orc2_unsigned, COMMUTE(orc1), ORC2,) */
842 /* DEFINE_BITWISE_UNSIGNED (bignum_bitwise_not1_unsigned, COMMUTE(not1), ARG2, bignum_bitwise_not) */
843 /* DEFINE_BITWISE_UNSIGNED (bignum_bitwise_orc1_unsigned, COMMUTE(orc2), ORC2,) */
844 DEFINE_BITWISE_UNSIGNED (bignum_bitwise_nand_unsigned, SWAP, AND, bignum_bitwise_not)
845 /* DEFINE_BITWISE_UNSIGNED (bignum_bitwise_set_unsigned, SWAP, CLEAR, bignum_bitwise_not) */
846
847 static bignum_type
848 bignum_bitwise_orc2_unsigned (bignum_type x, bignum_type y)
849 {
850 return (bignum_bitwise_not (bignum_bitwise_andc1 (x, y)));
851 }
852
853 static bignum_type
bignum_bitwise_orc1_unsigned(bignum_type x,bignum_type y)854 bignum_bitwise_orc1_unsigned (bignum_type x, bignum_type y)
855 {
856 return (bignum_bitwise_not (bignum_bitwise_andc2 (x, y)));
857 }
858
859 #undef SWAP
860 #undef COMMUTE
861 #undef DEFINE_BITWISE_UNSIGNED
862
863 #define DEFINE_BITWISE(NAME, X0, Y0, PP, PN, NP, NN) \
864 bignum_type \
865 NAME (bignum_type x, bignum_type y) \
866 { \
867 if (BIGNUM_ZERO_P (x)) return (Y0); \
868 if (BIGNUM_ZERO_P (y)) return (X0); \
869 return \
870 ((BIGNUM_POSITIVE_P (x)) \
871 ? ((BIGNUM_POSITIVE_P (y)) \
872 ? (bignum_bitwise_##PP##_unsigned (x, y)) \
873 : (bignum_bitwise_##PN##_unsigned (x, (bignum_bitwise_not (y))))) \
874 : ((BIGNUM_POSITIVE_P (y)) \
875 ? (bignum_bitwise_##NP##_unsigned ((bignum_bitwise_not (x)), y)) \
876 : (bignum_bitwise_##NN##_unsigned \
877 ((bignum_bitwise_not (x)), (bignum_bitwise_not (y)))))); \
878 }
879
880 DEFINE_BITWISE (bignum_bitwise_and, (BIGNUM_ZERO ()), (BIGNUM_ZERO ()),
881 and, andc2, andc1, nor)
882
883 DEFINE_BITWISE (bignum_bitwise_andc2, x, (BIGNUM_ZERO ()),
884 andc2, and, nor, andc1)
885
886 DEFINE_BITWISE (bignum_bitwise_andc1, (BIGNUM_ZERO ()), y,
887 andc1, nor, and, andc2)
888
DEFINE_BITWISE(bignum_bitwise_xor,x,y,xor,eqv,eqv,xor)889 DEFINE_BITWISE (bignum_bitwise_xor, x, y, xor, eqv, eqv, xor)
890
891 DEFINE_BITWISE (bignum_bitwise_ior, x, y, ior, orc2, orc1, nand)
892
893 DEFINE_BITWISE (bignum_bitwise_nor,
894 (bignum_bitwise_not (x)),
895 (bignum_bitwise_not (y)),
896 nor, andc1, andc2, and)
897
898 DEFINE_BITWISE (bignum_bitwise_eqv,
899 (bignum_bitwise_not (x)),
900 (bignum_bitwise_not (y)),
901 eqv, xor, xor, eqv)
902
903 DEFINE_BITWISE (bignum_bitwise_orc2,
904 (BIGNUM_ONE (1)), (bignum_bitwise_not (y)),
905 orc2, ior, nand, orc1)
906
907 DEFINE_BITWISE (bignum_bitwise_orc1,
908 (bignum_bitwise_not (x)), (BIGNUM_ONE (1)),
909 orc1, nand, ior, orc2)
910
911 DEFINE_BITWISE (bignum_bitwise_nand, (BIGNUM_ONE (1)), (BIGNUM_ONE (1)),
912 nand, orc1, orc2, ior)
913
914 #undef DEFINE_BITWISE
915
916 /* General bit-twiddlers. */
917
918 /* (edit a a-pos b b-pos selector size)
919 = (ior (and a (shift-left selector a-pos))
920 (shift-left (and (shift-right b b-pos) (not selector)) a-pos))
921
922 In other words, replace a[a-pos + i] by b[b-pos + i] if selector[i]
923 is zero, and shift the result right by pos (which can be negative).
924 For example, (extract size pos x) = (edit x pos 0 0 (mask size 0)). */
925
926 #if 0
927 bignum_type
928 bignum_edit_bit_field (bignum_type x, unsigned long x_position,
929 bignum_type y, unsigned long y_position,
930 bignum_type selector, unsigned long size)
931 {
932 BIGNUM_ASSERT (!BIGNUM_NEGATIVE_P (x));
933 BIGNUM_ASSERT (!BIGNUM_NEGATIVE_P (y));
934 BIGNUM_ASSERT (!BIGNUM_NEGATIVE_P (selector));
935 if (BIGNUM_ZERO_P (selector))
936 return (x);
937 {
938 bignum_length_type x_length = (BIGNUM_LENGTH (x));
939 bignum_length_type y_length = (BIGNUM_LENGTH (y));
940 bignum_length_type selector_length = (BIGNUM_LENGTH (selector));
941 bignum_length_type r_length = (MAX (x_length, (y_position + size)));
942 bignum_type r = (bignum_allocate (r_length, (BIGNUM_NEGATIVE_P (x))));
943 bignum_length_type x_digit_position = (x_position / BIGNUM_DIGIT_LENGTH);
944 bignum_length_type y_digit_position = (y_position / BIGNUM_DIGIT_LENGTH);
945 bignum_length_type x_bit_position = (x_position % BIGNUM_DIGIT_LENGTH);
946 bignum_length_type y_bit_position = (y_position % BIGNUM_DIGIT_LENGTH);
947 bignum_digit_type *x_scan = (BIGNUM_START_PTR (x));
948 bignum_digit_type *y_scan = (BIGNUM_START_PTR (y));
949 bignum_digit_type *selector_scan = (BIGNUM_START_PTR (selector));
950 bignum_digit_type *r_scan = (BIGNUM_START_PTR (r));
951 bignum_digit_type *x_end = (x_scan + x_length);
952 bignum_digit_type *y_end = (y_scan + y_length);
953 bignum_digit_type *selector_end = (selector_scan + selector_length);
954 {
955 bignum_digit_type *stop = (x_scan + (MIN (x_digit_position, x_length)));
956 while (x_scan < stop)
957 (*r_scan++) = (*x_scan++);
958 }
959 /* Four cases to deal with, depending on whether or not x and y
960 have more digits. */
961 y_scan += (MIN (y_digit_position, (y_end - y_scan)));
962 if ((x_scan < x_end) && (y_scan < y_end))
963 {
964 bignum_digit_type x_low_mask = (BIGNUM_DIGIT_ONES (x_bit_position));
965 bignum_digit_type x_high_mask
966 = (BIGNUM_DIGIT_ONES (BIGNUM_DIGIT_LENGTH - x_bit_position));
967 bignum_digit_type y_low_mask = (BIGNUM_DIGIT_ONES (y_bit_position));
968 bignum_digit_type y_high_mask
969 = (BIGNUM_DIGIT_ONES (BIGNUM_DIGIT_LENGTH - y_bit_position));
970 bignum_digit_type r_carry = ((*x_scan) & x_low_mask);
971 bignum_digit_type x_carry
972 = (((*x_scan++) >> x_bit_position) & x_high_mask);
973 bignum_digit_type y_carry
974 = (((*y_scan++) >> y_bit_position) & y_high_mask);
975 while ((x_scan < x_end) && (y_scan < y_end))
976 {
977 bignum_digit_type selector = (*selector_scan++);
978 bignum_digit_type x = (*x_scan++);
979 bignum_digit_type y = (*y_scan++);
980 (*r_scan++)
981 = (r_carry
982 | ((x_carry ) << x_bit_position)
983
984
985 | ((((x >> x_bit_position) & x_high_mask &~ selector)
986 | (high_y & selector))
987 << x_bit_position));
988 carry = ...;
989 }
990 }
991 else if (x_scan < x_end)
992 {
993 }
994 else if (y_scan < y_end)
995 {
996 }
997 else
998 {
999 }
1000 }
1001 }
1002 #endif /* 0 */
1003
1004 /* (splice a a-pos b b-pos size)
1005 = (edit a a-pos b b-pos (mask size) size)
1006
1007 Thus, e.g., (extract size pos x) = (splice x pos 0 0 size). */
1008
1009 #if 0
1010 bignum_type
1011 bignum_splice_bit_field (bignum_type x, unsigned long x_position,
1012 bignum_type y, unsigned long y_position,
1013 unsigned long size)
1014 {
1015 ...
1016 }
1017 #endif /* 0 */
1018
1019 /* Ones-complement length. */
1020
1021 unsigned long
1022 bignum_length_in_bits (bignum_type bignum)
1023 {
1024 if (BIGNUM_ZERO_P (bignum))
1025 return (0);
1026 {
1027 bignum_length_type index = ((BIGNUM_LENGTH (bignum)) - 1);
1028 return
1029 ((ulong_length_in_bits (BIGNUM_REF (bignum, index)))
1030 + (BIGNUM_DIGIT_LENGTH * index));
1031 }
1032 }
1033
1034 /* Two's-complement length. */
1035
1036 unsigned long
bignum_integer_length(bignum_type bignum)1037 bignum_integer_length (bignum_type bignum)
1038 {
1039 unsigned long length_in_bits = (bignum_length_in_bits (bignum));
1040 if ((BIGNUM_ZERO_P (bignum)) || (BIGNUM_POSITIVE_P (bignum)))
1041 return (length_in_bits);
1042 /* else if (BIGNUM_NEGATIVE_P (bignum)) */
1043 {
1044 /* We have to test whether it is a negative power of two. If so,
1045 we treat its length as one less than bignum_length_in_bits,
1046 because we are really measuring the length of the finite
1047 sequence of bits before the infinite sequence of zero bits (for
1048 nonnegative integers) or one bits (for negative integers) in the
1049 integer's general two's-complement representation. Thus,
1050 negative powers of two appear to have one fewer bit. */
1051 bignum_digit_type *scan = (BIGNUM_START_PTR (bignum));
1052 bignum_digit_type *end = (scan + (BIGNUM_LENGTH (bignum)) - 1);
1053 while (scan < end)
1054 if (0 != (*scan++))
1055 return (length_in_bits);
1056 return (length_in_bits - (0 == ((*end) & ((*end) - 1))));
1057 }
1058 }
1059
1060 long
bignum_first_set_bit(bignum_type bignum)1061 bignum_first_set_bit (bignum_type bignum)
1062 {
1063 if (BIGNUM_ZERO_P (bignum))
1064 return (-1);
1065 {
1066 bignum_digit_type *start = (BIGNUM_START_PTR (bignum));
1067 bignum_digit_type *scan = start;
1068 bignum_digit_type *end = (scan + (BIGNUM_LENGTH (bignum)));
1069 while (scan < end)
1070 {
1071 if ((*scan) != 0) break;
1072 scan += 1;
1073 }
1074 BIGNUM_ASSERT (scan < end);
1075 return
1076 ((ulong_bit_count (((*scan) ^ ((*scan) - 1)) >> 1))
1077 + ((scan - start) * BIGNUM_DIGIT_LENGTH));
1078 }
1079 }
1080
1081 static inline unsigned long
digits_bit_count(bignum_digit_type ** scan_loc,unsigned long count)1082 digits_bit_count (bignum_digit_type **scan_loc, unsigned long count)
1083 {
1084 unsigned long bit_count = 0;
1085 bignum_digit_type *end = ((*scan_loc) + count);
1086 while ((*scan_loc) < end)
1087 bit_count += (ulong_bit_count (* ((*scan_loc) ++)));
1088 return (bit_count);
1089 }
1090
1091 static inline unsigned long
digits_hamming_distance(bignum_digit_type ** x_scan_loc,bignum_digit_type ** y_scan_loc,unsigned long count)1092 digits_hamming_distance (bignum_digit_type **x_scan_loc,
1093 bignum_digit_type **y_scan_loc,
1094 unsigned long count)
1095 {
1096 unsigned long hamming_distance = 0;
1097 bignum_digit_type *end = ((*x_scan_loc) + count);
1098 while ((*x_scan_loc) < end)
1099 hamming_distance
1100 += (ulong_bit_count ((* ((*x_scan_loc) ++)) ^ (* ((*y_scan_loc) ++))));
1101 return (hamming_distance);
1102 }
1103
1104 /* Hamming distance: the easy case. */
1105
1106 static unsigned long
bignum_positive_hamming_distance(bignum_type x,bignum_type y)1107 bignum_positive_hamming_distance (bignum_type x, bignum_type y)
1108 {
1109 BIGNUM_ASSERT (!BIGNUM_ZERO_P (x));
1110 BIGNUM_ASSERT (!BIGNUM_ZERO_P (y));
1111 BIGNUM_ASSERT (BIGNUM_POSITIVE_P (x));
1112 BIGNUM_ASSERT (BIGNUM_POSITIVE_P (y));
1113 if ((BIGNUM_LENGTH (x)) < (BIGNUM_LENGTH (y)))
1114 {
1115 bignum_type t = x; x = y; y = t;
1116 }
1117 {
1118 bignum_digit_type *x_scan = (BIGNUM_START_PTR (x));
1119 bignum_digit_type *y_scan = (BIGNUM_START_PTR (y));
1120 unsigned long hamming_distance
1121 = (digits_hamming_distance
1122 ((&x_scan), (&y_scan), (BIGNUM_LENGTH (y))));
1123 hamming_distance
1124 += (digits_bit_count
1125 ((&x_scan), ((BIGNUM_LENGTH (x)) - (BIGNUM_LENGTH (y)))));
1126 return (hamming_distance);
1127 }
1128 }
1129
1130 /* Hamming distance: the hard case. */
1131
1132 #if 0
1133 /* Is this actually faster than (hamming-distance (not x) (not y))? */
1134
1135 #define MIN(A,B) (((A) < (B)) ? (A) : (B))
1136
1137 static unsigned long
1138 bignum_negative_hamming_distance (bignum_type x, bignum_type y)
1139 {
1140 BIGNUM_ASSERT (!BIGNUM_ZERO_P (x));
1141 BIGNUM_ASSERT (!BIGNUM_ZERO_P (y));
1142 BIGNUM_ASSERT (BIGNUM_NEGATIVE_P (x));
1143 BIGNUM_ASSERT (BIGNUM_NEGATIVE_P (y));
1144 {
1145 bignum_digit_type *x_scan = (BIGNUM_START_PTR (x));
1146 bignum_digit_type *y_scan = (BIGNUM_START_PTR (y));
1147 bignum_digit_type *x_end = (x_scan + (BIGNUM_LENGTH (x)));
1148 bignum_digit_type *y_end = (y_scan + (BIGNUM_LENGTH (y)));
1149 bignum_digit_type x_digit, y_digit;
1150 unsigned long hamming_distance;
1151 /* Find the position of the first nonzero digit of x or y, and
1152 maybe exchange x and y to guarantee that x's is nonzero. */
1153 while (1)
1154 {
1155 BIGNUM_ASSERT (x_scan < x_end);
1156 BIGNUM_ASSERT (y_scan < y_end);
1157 x_digit = (*x_scan++);
1158 y_digit = (*y_scan++);
1159 if (x_digit != 0) break;
1160 if (y_digit != 0)
1161 {
1162 bignum_digit_type *t;
1163 t = x_scan; x_scan = y_scan; y_scan = t;
1164 t = x_end; x_end = y_end; y_end = t;
1165 x_digit = y_digit;
1166 y_digit = 0;
1167 break;
1168 }
1169 }
1170 /* Convert the first nonzero digit of x and the corresponding digit
1171 (possibly zero) of y to two's-complement. */
1172 x_digit = (-x_digit);
1173 y_digit = (-y_digit);
1174 hamming_distance
1175 = (ulong_bit_count ((x_digit ^ y_digit) & BIGNUM_DIGIT_MASK));
1176 /* Skip over zeroes in y. */
1177 if (y_digit == 0)
1178 {
1179 bignum_digit_type *y_ptr = y_scan;
1180 bignum_length_type zeroes;
1181 do {
1182 BIGNUM_ASSERT (y_scan < y_end);
1183 y_digit = (*y_scan++);
1184 } while (y_digit == 0);
1185 /* If we any more zeroes, compute the Hamming distance of that
1186 segment as if all corresponding digits of x were ones. */
1187 zeroes = (y_scan - y_ptr - 1);
1188 hamming_distance += (zeroes * BIGNUM_DIGIT_LENGTH);
1189 /* Then subtract the amount by which this overestimated. */
1190 hamming_distance
1191 -= (digits_bit_count ((&x_scan), (MIN ((x_end - x_scan), zeroes))));
1192 /* Convert the first nonzero digit of y to two's-complement --
1193 we can subtract 1 because it is nonzero and (semantically,
1194 if not in the type system) unsigned, and hence positive. */
1195 hamming_distance
1196 += (ulong_bit_count
1197 ((y_digit - 1) ^ ((x_scan < x_end) ? (*x_scan++) : 0)));
1198 }
1199 /* Finally, scan over the overlapping parts of x and y and then the
1200 non-overlapping high part of whichever is longer. */
1201 hamming_distance
1202 += (digits_hamming_distance
1203 ((&x_scan), (&y_scan), (MIN ((x_end - x_scan), (y_end - y_scan)))));
1204 hamming_distance
1205 += ((x_scan < x_end)
1206 ? (digits_bit_count ((&x_scan), (x_end - x_scan)))
1207 : (digits_bit_count ((&y_scan), (y_end - y_scan))));
1208 return (hamming_distance);
1209 }
1210 }
1211 #endif /* 0 */
1212
1213 static unsigned long
bignum_bit_count_unsigned(bignum_type x)1214 bignum_bit_count_unsigned (bignum_type x)
1215 {
1216 bignum_digit_type *scan = (BIGNUM_START_PTR (x));
1217 return (digits_bit_count ((&scan), (BIGNUM_LENGTH (x))));
1218 }
1219
1220 static unsigned long
bignum_negative_hamming_distance(bignum_type x,bignum_type y)1221 bignum_negative_hamming_distance (bignum_type x, bignum_type y)
1222 {
1223 x = (bignum_bitwise_not (x));
1224 y = (bignum_bitwise_not (y));
1225 if (BIGNUM_ZERO_P (x)) return (bignum_bit_count_unsigned (y));
1226 if (BIGNUM_ZERO_P (y)) return (bignum_bit_count_unsigned (x));
1227 return (bignum_positive_hamming_distance (x, y));
1228 }
1229
1230 unsigned long
bignum_bit_count(bignum_type x)1231 bignum_bit_count (bignum_type x) /* a.k.a. Hamming weight, pop count */
1232 {
1233 if (BIGNUM_ZERO_P (x))
1234 return (0);
1235 if (BIGNUM_NEGATIVE_P (x))
1236 /* return (bignum_negative_hamming_distance (x, (BIGNUM_ONE (1)))); */
1237 return (bignum_bit_count_unsigned (bignum_bitwise_not (x)));
1238 else
1239 return (bignum_bit_count_unsigned (x));
1240 }
1241
1242 long
bignum_hamming_distance(bignum_type x,bignum_type y)1243 bignum_hamming_distance (bignum_type x, bignum_type y)
1244 {
1245 if (x == y)
1246 return (0);
1247 if (BIGNUM_ZERO_P (x))
1248 return ((BIGNUM_NEGATIVE_P (y)) ? (-1) : (bignum_bit_count_unsigned (y)));
1249 if (BIGNUM_ZERO_P (y))
1250 return ((BIGNUM_NEGATIVE_P (x)) ? (-1) : (bignum_bit_count_unsigned (x)));
1251 return
1252 ((BIGNUM_POSITIVE_P (x))
1253 ? ((BIGNUM_POSITIVE_P (y))
1254 ? (bignum_positive_hamming_distance (x, y))
1255 : (-1))
1256 : ((BIGNUM_POSITIVE_P (y))
1257 ? (-1)
1258 : (bignum_negative_hamming_distance (x, y))));
1259 }
1260
1261 #if 0
1262
1263 bignum_type
1264 bignum_nonnegative_one_bits (unsigned long size, unsigned long position)
1265 {
1266 if (size == 0)
1267 return (BIGNUM_ZERO ());
1268 {
1269 unsigned long total = (size + position);
1270 bignum_length_type total_digits = (total / BIGNUM_DIGIT_LENGTH);
1271 bignum_length_type total_bits = (total % BIGNUM_DIGIT_LENGTH);
1272 bignum_length_type zero_digits = (position / BIGNUM_DIGIT_LENGTH);
1273 bignum_length_type zero_bits = (position % BIGNUM_DIGIT_LENGTH);
1274 bignum_length_type one_digits = (size / BIGNUM_DIGIT_LENGTH);
1275 bignum_length_type one_bits = (size % BIGNUM_DIGIT_LENGTH);
1276 bignum_length_type first_nonzero_bits
1277 = (BIGNUM_DIGIT_LENGTH - zero_bits);
1278 bignum_length_type first_one_bits
1279 = ((one_bits < first_nonzero_bits) ? one_bits
1280 : (one_bits - first_nonzero_bits));
1281 bignum_length_type last_one_bits = (one_bits - first_one_bits);
1282 bignum_length_type length = (total_digits + (0 != total_bits));
1283 bignum_type r = (bignum_allocate (length, 0));
1284 bignum_digit_type *r_scan = (BIGNUM_START_PTR (r));
1285 bignum_digit_type *r_zero_end = (r_scan + zero_digits);
1286 bignum_digit_type *r_one_end
1287 = (r_zero_end + (first_one_bits != 0) + one_digits);
1288 BIGNUM_ASSERT ((r_one_end + (last_one_bits != 0)) == (r_scan + length));
1289 while (r_scan < r_zero_end)
1290 (*r_scan++) = 0;
1291 if (first_one_bits != 0)
1292 (*r_scan++) = ((BIGNUM_DIGIT_ONES (first_one_bits)) << zero_bits);
1293 while (r_scan < r_one_end)
1294 (*r_scan++) = BIGNUM_DIGIT_MASK;
1295 if (last_one_bits != 0)
1296 (*r_scan++) = (BIGNUM_DIGIT_ONES (last_one_bits));
1297 return (r);
1298 }
1299 }
1300
1301 #endif
1302
1303 bignum_type
bignum_nonnegative_one_bits(unsigned long size,unsigned long position)1304 bignum_nonnegative_one_bits (unsigned long size, unsigned long position)
1305 {
1306 return
1307 (bignum_shift_left
1308 ((bignum_bitwise_not (bignum_shift_left ((BIGNUM_ONE (1)), size))),
1309 position));
1310 }
1311
1312 bignum_type
bignum_negative_zero_bits(unsigned long n,unsigned long m)1313 bignum_negative_zero_bits (unsigned long n, unsigned long m)
1314 {
1315 return (bignum_bitwise_not (bignum_nonnegative_one_bits (n, m)));
1316 }
1317
1318 static bignum_type
bignum_shift_right_unsigned(bignum_type n,unsigned long digits,unsigned long bits)1319 bignum_shift_right_unsigned (bignum_type n,
1320 unsigned long digits,
1321 unsigned long bits)
1322 {
1323 bignum_length_type n_length = (BIGNUM_LENGTH (n));
1324 bignum_digit_type *n_start = (BIGNUM_START_PTR (n));
1325 bignum_digit_type *n_scan = (n_start + digits);
1326 bignum_digit_type *n_end = (n_start + n_length);
1327 bignum_length_type r_length
1328 = (n_length - digits
1329 - ((n_start < n_end) && (0 == ((n_end[-1]) >> bits))));
1330 bignum_type r = (bignum_allocate (r_length, 0));
1331 bignum_digit_type *r_scan = (BIGNUM_START_PTR (r));
1332 if (bits == 0)
1333 while (n_scan < n_end)
1334 (*r_scan++) = (*n_scan++);
1335 else
1336 {
1337 bignum_digit_type mask = (BIGNUM_DIGIT_ONES (bits));
1338 bignum_digit_type shift = (BIGNUM_DIGIT_LENGTH - bits);
1339 bignum_digit_type extra = ((*n_scan++) >> bits);
1340 while (n_scan < n_end)
1341 {
1342 bignum_digit_type digit = (*n_scan++);
1343 (*r_scan++) = (((digit & mask) << shift) | extra);
1344 extra = (digit >> bits);
1345 }
1346 if (extra != 0)
1347 (*r_scan++) = extra;
1348 BIGNUM_ASSERT (r_scan == ((BIGNUM_START_PTR (r)) + r_length));
1349 }
1350 return (r);
1351 }
1352
1353 bignum_type
bignum_shift_right(bignum_type n,unsigned long m)1354 bignum_shift_right (bignum_type n, unsigned long m)
1355 {
1356 unsigned long digits = (m / BIGNUM_DIGIT_LENGTH);
1357 unsigned long bits = (m % BIGNUM_DIGIT_LENGTH);
1358
1359 if (digits >= (BIGNUM_LENGTH (n)))
1360 return ((BIGNUM_NEGATIVE_P (n)) ? (BIGNUM_ONE (1)) : (BIGNUM_ZERO ()));
1361
1362 if (BIGNUM_NEGATIVE_P (n))
1363 return
1364 (bignum_bitwise_not
1365 (bignum_shift_right_unsigned ((bignum_bitwise_not (n)), digits, bits)));
1366 else
1367 return (bignum_shift_right_unsigned (n, digits, bits));
1368 }
1369
1370 bignum_type
bignum_shift_left(bignum_type n,unsigned long m)1371 bignum_shift_left (bignum_type n, unsigned long m)
1372 {
1373 unsigned long ln = (BIGNUM_LENGTH (n));
1374 unsigned long delta = 0;
1375 if ((m == 0) || (BIGNUM_ZERO_P (n)))
1376 return (n);
1377
1378 delta = (ulong_length_in_bits (BIGNUM_REF (n, (ln - 1))));
1379
1380 {
1381 unsigned long zeroes = (m / BIGNUM_DIGIT_LENGTH);
1382 unsigned long shift = (m % BIGNUM_DIGIT_LENGTH);
1383 unsigned long ln2
1384 = (((ln - 1) + ((delta + m) / BIGNUM_DIGIT_LENGTH))
1385 + (((delta + m) % BIGNUM_DIGIT_LENGTH) != 0));
1386 bignum_type result = (bignum_allocate (ln2, (BIGNUM_NEGATIVE_P (n))));
1387 bignum_digit_type * scan_n = (BIGNUM_START_PTR (n));
1388 bignum_digit_type * end_n = (scan_n + ln);
1389 bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
1390 while ((zeroes--) > 0)
1391 (*scan_result++) = 0;
1392 if (shift == 0)
1393 while (scan_n < end_n)
1394 (*scan_result++) = (*scan_n++);
1395 else
1396 {
1397 unsigned long temp = 0;
1398 while (scan_n < end_n)
1399 {
1400 bignum_digit_type digit = (*scan_n++);
1401 (*scan_result++) = (((digit << shift) & BIGNUM_DIGIT_MASK) | temp);
1402 temp = (digit >> (BIGNUM_DIGIT_LENGTH - shift));
1403 }
1404 if (temp != 0)
1405 (*scan_result) = temp;
1406 }
1407 return (result);
1408 }
1409 }
1410
1411 bignum_type
unsigned_long_to_shifted_bignum(unsigned long n,unsigned long m,int sign)1412 unsigned_long_to_shifted_bignum (unsigned long n, unsigned long m, int sign)
1413 {
1414 unsigned long delta = 0;
1415 if (n == 0)
1416 return (BIGNUM_ZERO ());
1417
1418 delta = (ulong_length_in_bits (n));
1419
1420 {
1421 unsigned long zeroes = (m / BIGNUM_DIGIT_LENGTH);
1422 unsigned long shift = (m % BIGNUM_DIGIT_LENGTH);
1423 unsigned long ln
1424 = (((delta + m) / BIGNUM_DIGIT_LENGTH)
1425 + (((delta + m) % BIGNUM_DIGIT_LENGTH) != 0));
1426 bignum_type result = (bignum_allocate (ln, sign));
1427 bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
1428 while ((zeroes--) > 0)
1429 (*scan_result++) = 0;
1430 (*scan_result++) = ((n << shift) & BIGNUM_DIGIT_MASK);
1431 n >>= (BIGNUM_DIGIT_LENGTH - shift);
1432 while (n > 0)
1433 {
1434 (*scan_result++) = (n & BIGNUM_DIGIT_MASK);
1435 n >>= BIGNUM_DIGIT_LENGTH;
1436 }
1437 return (result);
1438 }
1439 }
1440
1441 bignum_type
digit_stream_to_bignum(unsigned int n_digits,unsigned int (* producer)(bignum_procedure_context),bignum_procedure_context context,unsigned int radix,int negative_p)1442 digit_stream_to_bignum (unsigned int n_digits,
1443 unsigned int (*producer) (bignum_procedure_context),
1444 bignum_procedure_context context,
1445 unsigned int radix,
1446 int negative_p)
1447 {
1448 BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
1449 if (n_digits == 0)
1450 return (BIGNUM_ZERO ());
1451 if (n_digits == 1)
1452 {
1453 long digit = ((long) ((*producer) (context)));
1454 return (long_to_bignum (negative_p ? (- digit) : digit));
1455 }
1456 {
1457 /* This length will be at least as large as needed. */
1458 bignum_length_type length
1459 = (BIGNUM_BITS_TO_DIGITS
1460 (n_digits * (ulong_length_in_bits (radix))));
1461 {
1462 bignum_type result = (bignum_allocate_zeroed (length, negative_p));
1463 while ((n_digits--) > 0)
1464 {
1465 bignum_destructive_scale_up (result, ((bignum_digit_type) radix));
1466 bignum_destructive_add
1467 (result, ((bignum_digit_type) ((*producer) (context))));
1468 }
1469 return (bignum_trim (result));
1470 }
1471 }
1472 }
1473
1474 void
bignum_to_digit_stream(bignum_type bignum,unsigned int radix,void (* consumer)(bignum_procedure_context,long),bignum_procedure_context context)1475 bignum_to_digit_stream (bignum_type bignum,
1476 unsigned int radix,
1477 void (*consumer) (bignum_procedure_context, long),
1478 bignum_procedure_context context)
1479 {
1480 BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
1481 if (! (BIGNUM_ZERO_P (bignum)))
1482 {
1483 bignum_type working_copy = (bignum_copy (bignum));
1484 bignum_digit_type * start = (BIGNUM_START_PTR (working_copy));
1485 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (working_copy)));
1486 while (start < scan)
1487 {
1488 if ((scan[-1]) == 0)
1489 scan -= 1;
1490 else
1491 (*consumer)
1492 (context, (bignum_destructive_scale_down (working_copy, radix)));
1493 }
1494 BIGNUM_DEALLOCATE (working_copy);
1495 }
1496 return;
1497 }
1498
1499 long
bignum_max_digit_stream_radix(void)1500 bignum_max_digit_stream_radix (void)
1501 {
1502 return (BIGNUM_RADIX_ROOT);
1503 }
1504
1505 /* Comparisons */
1506
1507 static int
bignum_equal_p_unsigned(bignum_type x,bignum_type y)1508 bignum_equal_p_unsigned (bignum_type x, bignum_type y)
1509 {
1510 bignum_length_type length = (BIGNUM_LENGTH (x));
1511 if (length != ((bignum_length_type) (BIGNUM_LENGTH (y))))
1512 return (0);
1513 else
1514 {
1515 bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
1516 bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
1517 bignum_digit_type * end_x = (scan_x + length);
1518 while (scan_x < end_x)
1519 if ((*scan_x++) != (*scan_y++))
1520 return (0);
1521 return (1);
1522 }
1523 }
1524
1525 static enum bignum_comparison
bignum_compare_unsigned(bignum_type x,bignum_type y)1526 bignum_compare_unsigned (bignum_type x, bignum_type y)
1527 {
1528 bignum_length_type x_length = (BIGNUM_LENGTH (x));
1529 bignum_length_type y_length = (BIGNUM_LENGTH (y));
1530 if (x_length < y_length)
1531 return (bignum_comparison_less);
1532 if (x_length > y_length)
1533 return (bignum_comparison_greater);
1534 {
1535 bignum_digit_type * start_x = (BIGNUM_START_PTR (x));
1536 bignum_digit_type * scan_x = (start_x + x_length);
1537 bignum_digit_type * scan_y = ((BIGNUM_START_PTR (y)) + y_length);
1538 while (start_x < scan_x)
1539 {
1540 bignum_digit_type digit_x = (*--scan_x);
1541 bignum_digit_type digit_y = (*--scan_y);
1542 if (digit_x < digit_y)
1543 return (bignum_comparison_less);
1544 if (digit_x > digit_y)
1545 return (bignum_comparison_greater);
1546 }
1547 }
1548 return (bignum_comparison_equal);
1549 }
1550
1551 /* Addition */
1552
1553 static bignum_type
bignum_add_unsigned(bignum_type x,bignum_type y,int negative_p)1554 bignum_add_unsigned (bignum_type x, bignum_type y, int negative_p)
1555 {
1556 if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
1557 {
1558 bignum_type z = x;
1559 x = y;
1560 y = z;
1561 }
1562 {
1563 bignum_length_type x_length = (BIGNUM_LENGTH (x));
1564 bignum_type r = (bignum_allocate ((x_length + 1), negative_p));
1565 bignum_digit_type sum;
1566 bignum_digit_type carry = 0;
1567 bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
1568 bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
1569 {
1570 bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
1571 bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
1572 while (scan_y < end_y)
1573 {
1574 sum = ((*scan_x++) + (*scan_y++) + carry);
1575 if (sum < BIGNUM_RADIX)
1576 {
1577 (*scan_r++) = sum;
1578 carry = 0;
1579 }
1580 else
1581 {
1582 (*scan_r++) = (sum - BIGNUM_RADIX);
1583 carry = 1;
1584 }
1585 }
1586 }
1587 {
1588 bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
1589 if (carry != 0)
1590 while (scan_x < end_x)
1591 {
1592 sum = ((*scan_x++) + 1);
1593 if (sum < BIGNUM_RADIX)
1594 {
1595 (*scan_r++) = sum;
1596 carry = 0;
1597 break;
1598 }
1599 else
1600 (*scan_r++) = (sum - BIGNUM_RADIX);
1601 }
1602 while (scan_x < end_x)
1603 (*scan_r++) = (*scan_x++);
1604 }
1605 if (carry != 0)
1606 {
1607 (*scan_r) = 1;
1608 return (r);
1609 }
1610 return (bignum_shorten_length (r, x_length));
1611 }
1612 }
1613
1614 /* Subtraction */
1615
1616 static bignum_type
bignum_subtract_unsigned(bignum_type x,bignum_type y)1617 bignum_subtract_unsigned (bignum_type x, bignum_type y)
1618 {
1619 int negative_p = 0;
1620 switch (bignum_compare_unsigned (x, y))
1621 {
1622 case bignum_comparison_equal:
1623 return (BIGNUM_ZERO ());
1624 case bignum_comparison_less:
1625 {
1626 bignum_type z = x;
1627 x = y;
1628 y = z;
1629 }
1630 negative_p = 1;
1631 break;
1632 case bignum_comparison_greater:
1633 negative_p = 0;
1634 break;
1635 }
1636 {
1637 bignum_length_type x_length = (BIGNUM_LENGTH (x));
1638 bignum_type r = (bignum_allocate (x_length, negative_p));
1639 bignum_digit_type difference;
1640 bignum_digit_type borrow = 0;
1641 bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
1642 bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
1643 {
1644 bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
1645 bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
1646 while (scan_y < end_y)
1647 {
1648 difference = (((*scan_x++) - (*scan_y++)) - borrow);
1649 if (difference < 0)
1650 {
1651 (*scan_r++) = (difference + BIGNUM_RADIX);
1652 borrow = 1;
1653 }
1654 else
1655 {
1656 (*scan_r++) = difference;
1657 borrow = 0;
1658 }
1659 }
1660 }
1661 {
1662 bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
1663 if (borrow != 0)
1664 while (scan_x < end_x)
1665 {
1666 difference = ((*scan_x++) - borrow);
1667 if (difference < 0)
1668 (*scan_r++) = (difference + BIGNUM_RADIX);
1669 else
1670 {
1671 (*scan_r++) = difference;
1672 borrow = 0;
1673 break;
1674 }
1675 }
1676 BIGNUM_ASSERT (borrow == 0);
1677 while (scan_x < end_x)
1678 (*scan_r++) = (*scan_x++);
1679 }
1680 return (bignum_trim (r));
1681 }
1682 }
1683
1684 /* Multiplication
1685 Maximum value for product_low or product_high:
1686 ((R * R) + (R * (R - 2)) + (R - 1))
1687 Maximum value for carry: ((R * (R - 1)) + (R - 1))
1688 where R == BIGNUM_RADIX_ROOT */
1689
1690 static bignum_type
bignum_multiply_unsigned(bignum_type x,bignum_type y,int negative_p)1691 bignum_multiply_unsigned (bignum_type x, bignum_type y, int negative_p)
1692 {
1693 if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
1694 {
1695 bignum_type z = x;
1696 x = y;
1697 y = z;
1698 }
1699 {
1700 bignum_digit_type carry;
1701 bignum_digit_type y_digit_low;
1702 bignum_digit_type y_digit_high;
1703 bignum_digit_type x_digit_low;
1704 bignum_digit_type x_digit_high;
1705 bignum_digit_type product_low;
1706 bignum_digit_type * scan_r;
1707 bignum_digit_type * scan_y;
1708 bignum_length_type x_length = (BIGNUM_LENGTH (x));
1709 bignum_length_type y_length = (BIGNUM_LENGTH (y));
1710 bignum_type r =
1711 (bignum_allocate_zeroed ((x_length + y_length), negative_p));
1712 bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
1713 bignum_digit_type * end_x = (scan_x + x_length);
1714 bignum_digit_type * start_y = (BIGNUM_START_PTR (y));
1715 bignum_digit_type * end_y = (start_y + y_length);
1716 bignum_digit_type * start_r = (BIGNUM_START_PTR (r));
1717 #define x_digit x_digit_high
1718 #define y_digit y_digit_high
1719 #define product_high carry
1720 while (scan_x < end_x)
1721 {
1722 x_digit = (*scan_x++);
1723 x_digit_low = (HD_LOW (x_digit));
1724 x_digit_high = (HD_HIGH (x_digit));
1725 carry = 0;
1726 scan_y = start_y;
1727 scan_r = (start_r++);
1728 while (scan_y < end_y)
1729 {
1730 y_digit = (*scan_y++);
1731 y_digit_low = (HD_LOW (y_digit));
1732 y_digit_high = (HD_HIGH (y_digit));
1733 product_low =
1734 ((*scan_r) +
1735 (x_digit_low * y_digit_low) +
1736 (HD_LOW (carry)));
1737 product_high =
1738 ((x_digit_high * y_digit_low) +
1739 (x_digit_low * y_digit_high) +
1740 (HD_HIGH (product_low)) +
1741 (HD_HIGH (carry)));
1742 (*scan_r++) =
1743 (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
1744 carry =
1745 ((x_digit_high * y_digit_high) +
1746 (HD_HIGH (product_high)));
1747 }
1748 (*scan_r) += carry;
1749 }
1750 return (bignum_trim (r));
1751 #undef x_digit
1752 #undef y_digit
1753 #undef product_high
1754 }
1755 }
1756
1757 static bignum_type
bignum_multiply_unsigned_small_factor(bignum_type x,bignum_digit_type y,int negative_p)1758 bignum_multiply_unsigned_small_factor (bignum_type x, bignum_digit_type y,
1759 int negative_p)
1760 {
1761 bignum_length_type length_x = (BIGNUM_LENGTH (x));
1762 bignum_type p = (bignum_allocate ((length_x + 1), negative_p));
1763 bignum_destructive_copy (x, p);
1764 (BIGNUM_REF (p, length_x)) = 0;
1765 bignum_destructive_scale_up (p, y);
1766 return (bignum_trim (p));
1767 }
1768
1769 static void
bignum_destructive_scale_up(bignum_type bignum,bignum_digit_type factor)1770 bignum_destructive_scale_up (bignum_type bignum, bignum_digit_type factor)
1771 {
1772 bignum_digit_type carry = 0;
1773 bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
1774 bignum_digit_type two_digits;
1775 bignum_digit_type product_low;
1776 #define product_high carry
1777 bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum)));
1778 BIGNUM_ASSERT ((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
1779 while (scan < end)
1780 {
1781 two_digits = (*scan);
1782 product_low = ((factor * (HD_LOW (two_digits))) + (HD_LOW (carry)));
1783 product_high =
1784 ((factor * (HD_HIGH (two_digits))) +
1785 (HD_HIGH (product_low)) +
1786 (HD_HIGH (carry)));
1787 (*scan++) = (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
1788 carry = (HD_HIGH (product_high));
1789 }
1790 /* A carry here would be an overflow, i.e. it would not fit.
1791 Hopefully the callers allocate enough space that this will
1792 never happen.
1793 */
1794 BIGNUM_ASSERT (carry == 0);
1795 return;
1796 #undef product_high
1797 }
1798
1799 static void
bignum_destructive_add(bignum_type bignum,bignum_digit_type n)1800 bignum_destructive_add (bignum_type bignum, bignum_digit_type n)
1801 {
1802 bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
1803 bignum_digit_type digit;
1804 digit = ((*scan) + n);
1805 if (digit < BIGNUM_RADIX)
1806 {
1807 (*scan) = digit;
1808 return;
1809 }
1810 (*scan++) = (digit - BIGNUM_RADIX);
1811 while (1)
1812 {
1813 digit = ((*scan) + 1);
1814 if (digit < BIGNUM_RADIX)
1815 {
1816 (*scan) = digit;
1817 return;
1818 }
1819 (*scan++) = (digit - BIGNUM_RADIX);
1820 }
1821 }
1822
1823 /* Division */
1824
1825 /* For help understanding this algorithm, see:
1826 Knuth, Donald E., "The Art of Computer Programming",
1827 volume 2, "Seminumerical Algorithms"
1828 section 4.3.1, "Multiple-Precision Arithmetic". */
1829
1830 static void
bignum_divide_unsigned_large_denominator(bignum_type numerator,bignum_type denominator,bignum_type * quotient,bignum_type * remainder,int q_negative_p,int r_negative_p)1831 bignum_divide_unsigned_large_denominator (bignum_type numerator,
1832 bignum_type denominator,
1833 bignum_type * quotient,
1834 bignum_type * remainder,
1835 int q_negative_p,
1836 int r_negative_p)
1837 {
1838 bignum_length_type length_n = ((BIGNUM_LENGTH (numerator)) + 1);
1839 bignum_length_type length_d = (BIGNUM_LENGTH (denominator));
1840 bignum_type q =
1841 ((quotient != ((bignum_type *) 0))
1842 ? (bignum_allocate ((length_n - length_d), q_negative_p))
1843 : BIGNUM_OUT_OF_BAND);
1844 bignum_type u = (bignum_allocate (length_n, r_negative_p));
1845 int shift = 0;
1846 BIGNUM_ASSERT (length_d > 1);
1847 {
1848 bignum_digit_type v1 = (BIGNUM_REF ((denominator), (length_d - 1)));
1849 while (v1 < (BIGNUM_RADIX / 2))
1850 {
1851 v1 <<= 1;
1852 shift += 1;
1853 }
1854 }
1855 if (shift == 0)
1856 {
1857 bignum_destructive_copy (numerator, u);
1858 (BIGNUM_REF (u, (length_n - 1))) = 0;
1859 bignum_divide_unsigned_normalized (u, denominator, q);
1860 }
1861 else
1862 {
1863 bignum_type v = (bignum_allocate (length_d, 0));
1864 bignum_destructive_normalization (numerator, u, shift);
1865 bignum_destructive_normalization (denominator, v, shift);
1866 bignum_divide_unsigned_normalized (u, v, q);
1867 BIGNUM_DEALLOCATE (v);
1868 if (remainder != ((bignum_type *) 0))
1869 bignum_destructive_unnormalization (u, shift);
1870 }
1871 if (quotient != ((bignum_type *) 0))
1872 (*quotient) = (bignum_trim (q));
1873 if (remainder != ((bignum_type *) 0))
1874 (*remainder) = (bignum_trim (u));
1875 else
1876 BIGNUM_DEALLOCATE (u);
1877 return;
1878 }
1879
1880 static void
bignum_divide_unsigned_normalized(bignum_type u,bignum_type v,bignum_type q)1881 bignum_divide_unsigned_normalized (bignum_type u, bignum_type v, bignum_type q)
1882 {
1883 bignum_length_type u_length = (BIGNUM_LENGTH (u));
1884 bignum_length_type v_length = (BIGNUM_LENGTH (v));
1885 bignum_digit_type * u_start = (BIGNUM_START_PTR (u));
1886 bignum_digit_type * u_scan = (u_start + u_length);
1887 bignum_digit_type * u_scan_limit = (u_start + v_length);
1888 bignum_digit_type * u_scan_start = (u_scan - v_length);
1889 bignum_digit_type * v_start = (BIGNUM_START_PTR (v));
1890 bignum_digit_type * v_end = (v_start + v_length);
1891 bignum_digit_type * q_scan = 0;
1892 bignum_digit_type v1 = (v_end[-1]);
1893 bignum_digit_type v2 = (v_end[-2]);
1894 bignum_digit_type ph; /* high half of double-digit product */
1895 bignum_digit_type pl; /* low half of double-digit product */
1896 bignum_digit_type guess;
1897 bignum_digit_type gh; /* high half-digit of guess */
1898 bignum_digit_type ch; /* high half of double-digit comparand */
1899 bignum_digit_type v2l = (HD_LOW (v2));
1900 bignum_digit_type v2h = (HD_HIGH (v2));
1901 bignum_digit_type cl; /* low half of double-digit comparand */
1902 #define gl ph /* low half-digit of guess */
1903 #define uj pl
1904 #define qj ph
1905 bignum_digit_type gm; /* memory loc for reference parameter */
1906 if (q != BIGNUM_OUT_OF_BAND)
1907 q_scan = ((BIGNUM_START_PTR (q)) + (BIGNUM_LENGTH (q)));
1908 while (u_scan_limit < u_scan)
1909 {
1910 uj = (*--u_scan);
1911 if (uj != v1)
1912 {
1913 /* comparand =
1914 (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
1915 guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
1916 cl = (u_scan[-2]);
1917 ch = (bignum_digit_divide (uj, (u_scan[-1]), v1, (&gm)));
1918 guess = gm;
1919 }
1920 else
1921 {
1922 cl = (u_scan[-2]);
1923 ch = ((u_scan[-1]) + v1);
1924 guess = (BIGNUM_RADIX - 1);
1925 }
1926 while (1)
1927 {
1928 /* product = (guess * v2); */
1929 gl = (HD_LOW (guess));
1930 gh = (HD_HIGH (guess));
1931 pl = (v2l * gl);
1932 ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH (pl)));
1933 pl = (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl))));
1934 ph = ((v2h * gh) + (HD_HIGH (ph)));
1935 /* if (comparand >= product) */
1936 if ((ch > ph) || ((ch == ph) && (cl >= pl)))
1937 break;
1938 guess -= 1;
1939 /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
1940 ch += v1;
1941 /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
1942 if (ch >= BIGNUM_RADIX)
1943 break;
1944 }
1945 qj = (bignum_divide_subtract (v_start, v_end, guess, (--u_scan_start)));
1946 if (q != BIGNUM_OUT_OF_BAND)
1947 (*--q_scan) = qj;
1948 }
1949 return;
1950 #undef gl
1951 #undef uj
1952 #undef qj
1953 }
1954
1955 static bignum_digit_type
bignum_divide_subtract(bignum_digit_type * v_start,bignum_digit_type * v_end,bignum_digit_type guess,bignum_digit_type * u_start)1956 bignum_divide_subtract (bignum_digit_type * v_start,
1957 bignum_digit_type * v_end,
1958 bignum_digit_type guess,
1959 bignum_digit_type * u_start)
1960 {
1961 bignum_digit_type * v_scan = v_start;
1962 bignum_digit_type * u_scan = u_start;
1963 bignum_digit_type carry = 0;
1964 if (guess == 0) return (0);
1965 {
1966 bignum_digit_type gl = (HD_LOW (guess));
1967 bignum_digit_type gh = (HD_HIGH (guess));
1968 bignum_digit_type v;
1969 bignum_digit_type pl;
1970 bignum_digit_type vl;
1971 #define vh v
1972 #define ph carry
1973 #define diff pl
1974 while (v_scan < v_end)
1975 {
1976 v = (*v_scan++);
1977 vl = (HD_LOW (v));
1978 vh = (HD_HIGH (v));
1979 pl = ((vl * gl) + (HD_LOW (carry)));
1980 ph = ((vl * gh) + (vh * gl) + (HD_HIGH (pl)) + (HD_HIGH (carry)));
1981 diff = ((*u_scan) - (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl)))));
1982 if (diff < 0)
1983 {
1984 (*u_scan++) = (diff + BIGNUM_RADIX);
1985 carry = ((vh * gh) + (HD_HIGH (ph)) + 1);
1986 }
1987 else
1988 {
1989 (*u_scan++) = diff;
1990 carry = ((vh * gh) + (HD_HIGH (ph)));
1991 }
1992 }
1993 if (carry == 0)
1994 return (guess);
1995 diff = ((*u_scan) - carry);
1996 if (diff < 0)
1997 (*u_scan) = (diff + BIGNUM_RADIX);
1998 else
1999 {
2000 (*u_scan) = diff;
2001 return (guess);
2002 }
2003 #undef vh
2004 #undef ph
2005 #undef diff
2006 }
2007 /* Subtraction generated carry, implying guess is one too large.
2008 Add v back in to bring it back down. */
2009 v_scan = v_start;
2010 u_scan = u_start;
2011 carry = 0;
2012 while (v_scan < v_end)
2013 {
2014 bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
2015 if (sum < BIGNUM_RADIX)
2016 {
2017 (*u_scan++) = sum;
2018 carry = 0;
2019 }
2020 else
2021 {
2022 (*u_scan++) = (sum - BIGNUM_RADIX);
2023 carry = 1;
2024 }
2025 }
2026 if (carry == 1)
2027 {
2028 bignum_digit_type sum = ((*u_scan) + carry);
2029 (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
2030 }
2031 return (guess - 1);
2032 }
2033
2034 static void
bignum_divide_unsigned_medium_denominator(bignum_type numerator,bignum_digit_type denominator,bignum_type * quotient,bignum_type * remainder,int q_negative_p,int r_negative_p)2035 bignum_divide_unsigned_medium_denominator (bignum_type numerator,
2036 bignum_digit_type denominator,
2037 bignum_type * quotient,
2038 bignum_type * remainder,
2039 int q_negative_p,
2040 int r_negative_p)
2041 {
2042 bignum_length_type length_n = (BIGNUM_LENGTH (numerator));
2043 bignum_length_type length_q;
2044 bignum_type q;
2045 int shift = 0;
2046 /* Because `bignum_digit_divide' requires a normalized denominator. */
2047 while (denominator < (BIGNUM_RADIX / 2))
2048 {
2049 denominator <<= 1;
2050 shift += 1;
2051 }
2052 if (shift == 0)
2053 {
2054 length_q = length_n;
2055 q = (bignum_allocate (length_q, q_negative_p));
2056 bignum_destructive_copy (numerator, q);
2057 }
2058 else
2059 {
2060 length_q = (length_n + 1);
2061 q = (bignum_allocate (length_q, q_negative_p));
2062 bignum_destructive_normalization (numerator, q, shift);
2063 }
2064 {
2065 bignum_digit_type r = 0;
2066 bignum_digit_type * start = (BIGNUM_START_PTR (q));
2067 bignum_digit_type * scan = (start + length_q);
2068 bignum_digit_type qj;
2069 if (quotient != ((bignum_type *) 0))
2070 {
2071 while (start < scan)
2072 {
2073 r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
2074 (*scan) = qj;
2075 }
2076 (*quotient) = (bignum_trim (q));
2077 }
2078 else
2079 {
2080 while (start < scan)
2081 r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
2082 BIGNUM_DEALLOCATE (q);
2083 }
2084 if (remainder != ((bignum_type *) 0))
2085 {
2086 if (shift != 0)
2087 r >>= shift;
2088 (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
2089 }
2090 }
2091 return;
2092 }
2093
2094 static void
bignum_destructive_normalization(bignum_type source,bignum_type target,int shift_left)2095 bignum_destructive_normalization (bignum_type source, bignum_type target,
2096 int shift_left)
2097 {
2098 bignum_digit_type digit;
2099 bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
2100 bignum_digit_type carry = 0;
2101 bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
2102 bignum_digit_type * end_source = (scan_source + (BIGNUM_LENGTH (source)));
2103 bignum_digit_type * end_target = (scan_target + (BIGNUM_LENGTH (target)));
2104 int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
2105 bignum_digit_type mask = ((1UL << shift_right) - 1UL);
2106 while (scan_source < end_source)
2107 {
2108 digit = (*scan_source++);
2109 (*scan_target++) = (((digit & mask) << shift_left) | carry);
2110 carry = (digit >> shift_right);
2111 }
2112 if (scan_target < end_target)
2113 (*scan_target) = carry;
2114 else
2115 BIGNUM_ASSERT (carry == 0);
2116 return;
2117 }
2118
2119 static void
bignum_destructive_unnormalization(bignum_type bignum,int shift_right)2120 bignum_destructive_unnormalization (bignum_type bignum, int shift_right)
2121 {
2122 bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
2123 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
2124 bignum_digit_type digit;
2125 bignum_digit_type carry = 0;
2126 int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
2127 bignum_digit_type mask = ((1UL << shift_right) - 1UL);
2128 while (start < scan)
2129 {
2130 digit = (*--scan);
2131 (*scan) = ((digit >> shift_right) | carry);
2132 carry = ((digit & mask) << shift_left);
2133 }
2134 BIGNUM_ASSERT (carry == 0);
2135 }
2136
2137 /* This is a reduced version of the division algorithm, applied to the
2138 case of dividing two bignum digits by one bignum digit. It is
2139 assumed that the numerator and denominator are normalized. */
2140
2141 #define BDD_STEP(qn, j) \
2142 { \
2143 uj = (u[j]); \
2144 if (uj != v1) \
2145 { \
2146 uj_uj1 = (HD_CONS (uj, (u[j + 1]))); \
2147 guess = (uj_uj1 / v1); \
2148 comparand = (HD_CONS ((uj_uj1 % v1), (u[j + 2]))); \
2149 } \
2150 else \
2151 { \
2152 guess = (BIGNUM_RADIX_ROOT - 1); \
2153 comparand = (HD_CONS (((u[j + 1]) + v1), (u[j + 2]))); \
2154 } \
2155 while ((guess * v2) > comparand) \
2156 { \
2157 guess -= 1; \
2158 comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH); \
2159 if (comparand >= BIGNUM_RADIX) \
2160 break; \
2161 } \
2162 qn = (bignum_digit_divide_subtract (v1, v2, guess, (&u[j]))); \
2163 }
2164
2165 static bignum_digit_type
bignum_digit_divide(bignum_digit_type uh,bignum_digit_type ul,bignum_digit_type v,bignum_digit_type * q)2166 bignum_digit_divide (bignum_digit_type uh, bignum_digit_type ul,
2167 bignum_digit_type v, bignum_digit_type * q)
2168 {
2169 bignum_digit_type guess;
2170 bignum_digit_type comparand;
2171 bignum_digit_type v1;
2172 bignum_digit_type v2;
2173 bignum_digit_type uj;
2174 bignum_digit_type uj_uj1;
2175 bignum_digit_type q1;
2176 bignum_digit_type q2;
2177 bignum_digit_type u [4];
2178 if (uh == 0)
2179 {
2180 if (ul < v)
2181 {
2182 (*q) = 0;
2183 return (ul);
2184 }
2185 else if (ul == v)
2186 {
2187 (*q) = 1;
2188 return (0);
2189 }
2190 }
2191 (u[0]) = (HD_HIGH (uh));
2192 (u[1]) = (HD_LOW (uh));
2193 (u[2]) = (HD_HIGH (ul));
2194 (u[3]) = (HD_LOW (ul));
2195 v1 = (HD_HIGH (v));
2196 v2 = (HD_LOW (v));
2197 BDD_STEP (q1, 0);
2198 BDD_STEP (q2, 1);
2199 (*q) = (HD_CONS (q1, q2));
2200 return (HD_CONS ((u[2]), (u[3])));
2201 }
2202
2203 #undef BDD_STEP
2204
2205 #define BDDS_MULSUB(vn, un, carry_in) \
2206 { \
2207 product = ((vn * guess) + carry_in); \
2208 diff = (un - (HD_LOW (product))); \
2209 if (diff < 0) \
2210 { \
2211 un = (diff + BIGNUM_RADIX_ROOT); \
2212 carry = ((HD_HIGH (product)) + 1); \
2213 } \
2214 else \
2215 { \
2216 un = diff; \
2217 carry = (HD_HIGH (product)); \
2218 } \
2219 }
2220
2221 #define BDDS_ADD(vn, un, carry_in) \
2222 { \
2223 sum = (vn + un + carry_in); \
2224 if (sum < BIGNUM_RADIX_ROOT) \
2225 { \
2226 un = sum; \
2227 carry = 0; \
2228 } \
2229 else \
2230 { \
2231 un = (sum - BIGNUM_RADIX_ROOT); \
2232 carry = 1; \
2233 } \
2234 }
2235
2236 static bignum_digit_type
bignum_digit_divide_subtract(bignum_digit_type v1,bignum_digit_type v2,bignum_digit_type guess,bignum_digit_type * u)2237 bignum_digit_divide_subtract (bignum_digit_type v1, bignum_digit_type v2,
2238 bignum_digit_type guess, bignum_digit_type * u)
2239 {
2240 {
2241 bignum_digit_type product;
2242 bignum_digit_type diff;
2243 bignum_digit_type carry;
2244 BDDS_MULSUB (v2, (u[2]), 0);
2245 BDDS_MULSUB (v1, (u[1]), carry);
2246 if (carry == 0)
2247 return (guess);
2248 diff = ((u[0]) - carry);
2249 if (diff < 0)
2250 (u[0]) = (diff + BIGNUM_RADIX);
2251 else
2252 {
2253 (u[0]) = diff;
2254 return (guess);
2255 }
2256 }
2257 {
2258 bignum_digit_type sum;
2259 bignum_digit_type carry;
2260 BDDS_ADD(v2, (u[2]), 0);
2261 BDDS_ADD(v1, (u[1]), carry);
2262 if (carry == 1)
2263 (u[0]) += 1;
2264 }
2265 return (guess - 1);
2266 }
2267
2268 #undef BDDS_MULSUB
2269 #undef BDDS_ADD
2270
2271 static void
bignum_divide_unsigned_small_denominator(bignum_type numerator,bignum_digit_type denominator,bignum_type * quotient,bignum_type * remainder,int q_negative_p,int r_negative_p)2272 bignum_divide_unsigned_small_denominator (bignum_type numerator,
2273 bignum_digit_type denominator,
2274 bignum_type * quotient,
2275 bignum_type * remainder,
2276 int q_negative_p,
2277 int r_negative_p)
2278 {
2279 bignum_type q = (bignum_new_sign (numerator, q_negative_p));
2280 bignum_digit_type r = (bignum_destructive_scale_down (q, denominator));
2281 (*quotient) = (bignum_trim (q));
2282 if (remainder != ((bignum_type *) 0))
2283 (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
2284 return;
2285 }
2286
2287 /* Given (denominator > 1), it is fairly easy to show that
2288 (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
2289 that all digits are < BIGNUM_RADIX. */
2290
2291 static bignum_digit_type
bignum_destructive_scale_down(bignum_type bignum,bignum_digit_type denominator)2292 bignum_destructive_scale_down (bignum_type bignum,
2293 bignum_digit_type denominator)
2294 {
2295 bignum_digit_type numerator;
2296 bignum_digit_type remainder = 0;
2297 bignum_digit_type two_digits;
2298 #define quotient_high remainder
2299 bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
2300 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
2301 BIGNUM_ASSERT ((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
2302 while (start < scan)
2303 {
2304 two_digits = (*--scan);
2305 numerator = (HD_CONS (remainder, (HD_HIGH (two_digits))));
2306 quotient_high = (numerator / denominator);
2307 numerator = (HD_CONS ((numerator % denominator), (HD_LOW (two_digits))));
2308 (*scan) = (HD_CONS (quotient_high, (numerator / denominator)));
2309 remainder = (numerator % denominator);
2310 }
2311 return (remainder);
2312 #undef quotient_high
2313 }
2314
2315 static bignum_type
bignum_remainder_unsigned_small_denominator(bignum_type n,bignum_digit_type d,int negative_p)2316 bignum_remainder_unsigned_small_denominator (bignum_type n,
2317 bignum_digit_type d,
2318 int negative_p)
2319 {
2320 bignum_digit_type two_digits;
2321 bignum_digit_type * start = (BIGNUM_START_PTR (n));
2322 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (n)));
2323 bignum_digit_type r = 0;
2324 BIGNUM_ASSERT ((d > 1) && (d < BIGNUM_RADIX_ROOT));
2325 while (start < scan)
2326 {
2327 two_digits = (*--scan);
2328 r =
2329 ((HD_CONS (((HD_CONS (r, (HD_HIGH (two_digits)))) % d),
2330 (HD_LOW (two_digits))))
2331 % d);
2332 }
2333 return (bignum_digit_to_bignum (r, negative_p));
2334 }
2335
2336 static bignum_type
bignum_digit_to_bignum(bignum_digit_type digit,int negative_p)2337 bignum_digit_to_bignum (bignum_digit_type digit, int negative_p)
2338 {
2339 if (digit == 0)
2340 return (BIGNUM_ZERO ());
2341 else
2342 {
2343 bignum_type result = (bignum_allocate (1, negative_p));
2344 (BIGNUM_REF (result, 0)) = digit;
2345 return (result);
2346 }
2347 }
2348
2349 /* Allocation */
2350
2351 static bignum_type
bignum_allocate(bignum_length_type length,int negative_p)2352 bignum_allocate (bignum_length_type length, int negative_p)
2353 {
2354 BIGNUM_ASSERT ((0 <= length) && (length <= BIGNUM_LENGTH_MAX));
2355 {
2356 bignum_type result = (BIGNUM_ALLOCATE (length));
2357 BIGNUM_SET_HEADER (result, length, negative_p);
2358 return (result);
2359 }
2360 }
2361
2362 static bignum_type
bignum_allocate_zeroed(bignum_length_type length,int negative_p)2363 bignum_allocate_zeroed (bignum_length_type length, int negative_p)
2364 {
2365 BIGNUM_ASSERT ((0 <= length) && (length <= BIGNUM_LENGTH_MAX));
2366 {
2367 bignum_type result = (BIGNUM_ALLOCATE (length));
2368 bignum_digit_type * scan = (BIGNUM_START_PTR (result));
2369 bignum_digit_type * end = (scan + length);
2370 BIGNUM_SET_HEADER (result, length, negative_p);
2371 while (scan < end)
2372 (*scan++) = 0;
2373 return (result);
2374 }
2375 }
2376
2377 static bignum_type
bignum_shorten_length(bignum_type bignum,bignum_length_type length)2378 bignum_shorten_length (bignum_type bignum, bignum_length_type length)
2379 {
2380 bignum_length_type current_length = (BIGNUM_LENGTH (bignum));
2381 BIGNUM_ASSERT ((0 <= length) && (length <= current_length));
2382 if (length < current_length)
2383 {
2384 BIGNUM_SET_HEADER
2385 (bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
2386 BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
2387 }
2388 return (bignum);
2389 }
2390
2391 static bignum_type
bignum_trim(bignum_type bignum)2392 bignum_trim (bignum_type bignum)
2393 {
2394 bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
2395 bignum_digit_type * end = (start + (BIGNUM_LENGTH (bignum)));
2396 bignum_digit_type * scan = end;
2397 while ((start <= scan) && ((*--scan) == 0))
2398 ;
2399 scan += 1;
2400 if (scan < end)
2401 {
2402 bignum_length_type length = (scan - start);
2403 BIGNUM_SET_HEADER
2404 (bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
2405 BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
2406 }
2407 return (bignum);
2408 }
2409
2410 /* Copying */
2411
2412 static bignum_type
bignum_copy(bignum_type source)2413 bignum_copy (bignum_type source)
2414 {
2415 bignum_type target =
2416 (bignum_allocate ((BIGNUM_LENGTH (source)), (BIGNUM_NEGATIVE_P (source))));
2417 bignum_destructive_copy (source, target);
2418 return (target);
2419 }
2420
2421 static bignum_type
bignum_new_sign(bignum_type bignum,int negative_p)2422 bignum_new_sign (bignum_type bignum, int negative_p)
2423 {
2424 bignum_type result =
2425 (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
2426 bignum_destructive_copy (bignum, result);
2427 return (result);
2428 }
2429
2430 static bignum_type
bignum_maybe_new_sign(bignum_type bignum,int negative_p)2431 bignum_maybe_new_sign (bignum_type bignum, int negative_p)
2432 {
2433 #ifndef BIGNUM_FORCE_NEW_RESULTS
2434 if ((BIGNUM_NEGATIVE_P (bignum)) ? negative_p : (! negative_p))
2435 return (bignum);
2436 else
2437 #endif
2438 {
2439 bignum_type result =
2440 (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
2441 bignum_destructive_copy (bignum, result);
2442 return (result);
2443 }
2444 }
2445
2446 static void
bignum_destructive_copy(bignum_type source,bignum_type target)2447 bignum_destructive_copy (bignum_type source, bignum_type target)
2448 {
2449 bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
2450 bignum_digit_type * end_source =
2451 (scan_source + (BIGNUM_LENGTH (source)));
2452 bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
2453 while (scan_source < end_source)
2454 (*scan_target++) = (*scan_source++);
2455 return;
2456 }
2457