| blob | f96890aa5430332550c5af9c2542e82502a3486f |
1 ;;;; code to implement bignum support
3 ;;;; This software is part of the SBCL system. See the README file for
4 ;;;; more information.
5 ;;;;
6 ;;;; This software is derived from the CMU CL system, which was
7 ;;;; written at Carnegie Mellon University and released into the
8 ;;;; public domain. The software is in the public domain and is
9 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
10 ;;;; files for more information.
13 \f
14 ;;;; notes
16 ;;; comments from CMU CL:
17 ;;; These symbols define the interface to the number code:
18 ;;; add-bignums multiply-bignums negate-bignum subtract-bignum
19 ;;; multiply-bignum-and-fixnum multiply-fixnums
20 ;;; bignum-ashift-right bignum-ashift-left bignum-gcd
21 ;;; bignum-to-float bignum-integer-length
22 ;;; bignum-logical-and bignum-logical-ior bignum-logical-xor
23 ;;; bignum-logical-not bignum-load-byte bignum-deposit-byte
24 ;;; bignum-truncate bignum-plus-p bignum-compare make-small-bignum
25 ;;; bignum-logbitp bignum-logcount
26 ;;; These symbols define the interface to the compiler:
27 ;;; bignum-type bignum-element-type bignum-index %allocate-bignum
28 ;;; %bignum-length %bignum-set-length %bignum-ref %bignum-set
29 ;;; %digit-0-or-plusp %add-with-carry %subtract-with-borrow
30 ;;; %multiply-and-add %multiply %lognot %logand %logior %logxor
31 ;;; %fixnum-to-digit %floor %fixnum-digit-with-correct-sign %ashl
32 ;;; %ashr %digit-logical-shift-right))
34 ;;; The following interfaces will either be assembler routines or code
35 ;;; sequences expanded into the code as basic bignum operations:
36 ;;; General:
37 ;;; %BIGNUM-LENGTH
38 ;;; %ALLOCATE-BIGNUM
39 ;;; %BIGNUM-REF
40 ;;; %NORMALIZE-BIGNUM
41 ;;; %BIGNUM-SET-LENGTH
42 ;;; %FIXNUM-DIGIT-WITH-CORRECT-SIGN
43 ;;; %SIGN-DIGIT
44 ;;; %ASHR
45 ;;; %ASHL
46 ;;; %BIGNUM-0-OR-PLUSP
47 ;;; %DIGIT-LOGICAL-SHIFT-RIGHT
48 ;;; General (May not exist when done due to sole use in %-routines.)
49 ;;; %DIGIT-0-OR-PLUSP
50 ;;; Addition:
51 ;;; %ADD-WITH-CARRY
52 ;;; Subtraction:
53 ;;; %SUBTRACT-WITH-BORROW
54 ;;; Multiplication
55 ;;; %MULTIPLY
56 ;;; Negation
57 ;;; %LOGNOT
58 ;;; Shifting (in place)
59 ;;; %NORMALIZE-BIGNUM-BUFFER
60 ;;; GCD/Relational operators:
61 ;;; %DIGIT-COMPARE
62 ;;; %DIGIT-GREATER
63 ;;; Relational operators:
64 ;;; %LOGAND
65 ;;; %LOGIOR
66 ;;; %LOGXOR
67 ;;; LDB
68 ;;; %FIXNUM-TO-DIGIT
69 ;;; TRUNCATE
70 ;;; %FLOOR
71 ;;;
72 ;;; Note: The floating routines know about the float representation.
73 ;;;
74 ;;; PROBLEM 1:
75 ;;; There might be a problem with various LET's and parameters that take a
76 ;;; digit value. We need to write these so those things stay in machine
77 ;;; registers and number stack slots. I bind locals to these values, and I
78 ;;; use function on them -- ZEROP, ASH, etc.
79 ;;;
80 ;;; PROBLEM 2:
81 ;;; In shifting and byte operations, I use masks and logical operations that
82 ;;; could result in intermediate bignums. This is hidden by the current system,
83 ;;; but I may need to write these in a way that keeps these masks and logical
84 ;;; operations from diving into the Lisp level bignum code.
85 ;;;
86 ;;; To do:
87 ;;; fixnums
88 ;;; logior, logxor, logand
89 ;;; depending on relationals, < (twice) and <= (twice)
90 ;;; or write compare thing (twice).
91 ;;; LDB on fixnum with bignum result.
92 ;;; DPB on fixnum with bignum result.
93 ;;; TRUNCATE returns zero or one as one value and fixnum or minus fixnum
94 ;;; for the other value when given (truncate fixnum bignum).
95 ;;; Returns (truncate bignum fixnum) otherwise.
96 ;;; addition
97 ;;; subtraction (twice)
98 ;;; multiply
99 ;;; GCD
100 ;;; Write MASK-FIELD and DEPOSIT-FIELD in terms of logical operations.
101 ;;; DIVIDE
102 ;;; IF (/ x y) with bignums:
103 ;;; do the truncate, and if rem is 0, return quotient.
104 ;;; if rem is non-0
105 ;;; gcd of x and y.
106 ;;; "truncate" each by gcd, ignoring remainder 0.
107 ;;; form ratio of each result, bottom is positive.
108 \f
109 ;;;; What's a bignum?
117 \f
118 ;;;; internal inline routines
120 ;;; %ALLOCATE-BIGNUM must zero all elements.
125 ;;; Extract the length of the bignum.
130 ;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
131 ;;; to be able to return the digit somewhere no one looks for real objects.
142 ;;; Return T if digit is positive, or NIL if negative.
153 ;;; This should be in assembler, and should not cons intermediate
154 ;;; results. It returns a bignum digit and a carry resulting from adding
155 ;;; together a, b, and an incoming carry.
161 ;;; This should be in assembler, and should not cons intermediate
162 ;;; results. It returns a bignum digit and a borrow resulting from
163 ;;; subtracting b from a, and subtracting a possible incoming borrow.
164 ;;;
165 ;;; We really do: a - b - 1 + borrow, where borrow is either 0 or 1.
171 ;;; Multiply two digit-size numbers, returning a 2*digit-size result
172 ;;; split into two digit-size quantities.
177 ;;; This multiplies x-digit and y-digit, producing high and low digits
178 ;;; manifesting the result. Then it adds the low digit, res-digit, and
179 ;;; carry-in-digit. Any carries (note, you still have to add two digits
180 ;;; at a time possibly producing two carries) from adding these three
181 ;;; digits get added to the high digit from the multiply, producing the
182 ;;; next carry digit. Res-digit is optional since two uses of this
183 ;;; primitive multiplies a single digit bignum by a multiple digit
184 ;;; bignum, and in this situation there is no need for a result buffer
185 ;;; accumulating partial results which is where the res-digit comes
186 ;;; from.
196 ;;; Each of these does the digit-size unsigned op.
208 ;;; This takes a fixnum and sets it up as an unsigned digit-size
209 ;;; quantity.
215 ;;; This takes three digits and returns the FLOOR'ed result of
216 ;;; dividing the first two as a 2*digit-size integer by the third.
217 ;;;
218 ;;; Do weird LET and SETQ stuff to bamboozle the compiler into allowing
219 ;;; the %FLOOR transform to expand into pseudo-assembler for which the
220 ;;; compiler can later correctly allocate registers.
227 ;;; Convert the digit to a regular integer assuming that the digit is signed.
232 digit))
234 ;;; Do an arithmetic shift right of data even though bignum-element-type is
235 ;;; unsigned.
241 ;;; This takes a digit-size quantity and shifts it to the left,
242 ;;; returning a digit-size quantity.
248 ;;; Do an unsigned (logical) right shift of a digit by Count.
254 ;;; Change the length of bignum to be newlen. Newlen must be the same or
255 ;;; smaller than the old length, and any elements beyond newlen must be zeroed.
261 ;;; This returns 0 or "-1" depending on whether the bignum is positive. This
262 ;;; is suitable for infinite sign extension to complete additions,
263 ;;; subtractions, negations, etc. This cannot return a -1 represented as
264 ;;; a negative fixnum since it would then have to low zeros.
271 ;;; These take two digit-size quantities and compare or contrast them
272 ;;; without wasting time with incorrect type checking.
278 \f
280 \f
281 ;;;; addition
313 carry)))
316 ;;; This takes the longer of two bignums and propagates the carry through its
317 ;;; remaining high order digits.
333 \f
334 ;;;; subtraction
338 ;;; This subtracts b from a plugging result into res. Return-fun is the
339 ;;; function to call that fixes up the result returning any useful values, such
340 ;;; as the result. This macro may evaluate its arguments more than once.
376 ;;; Operations requiring a subtraction without the overhead of intermediate
377 ;;; results, such as GCD, use this. It assumes Result is big enough for the
378 ;;; result.
383 %normalize-bignum-buffer))
389 %normalize-bignum-buffer))
390 \f
391 ;;;; multiplication
414 (%multiply-and-add x
417 carry-digit)
463 low
470 \f
471 ;;;; BIGNUM-REPLACE and WITH-BIGNUM-BUFFERS
476 src
477 &key
479 end1
481 end2
482 from-end)
519 #!+sb-doc
520 "WITH-BIGNUM-BUFFERS ({(var size [init])}*) Form*"
535 \f
536 ;;;; GCD
539 ;; The asserts in the GCD implementation are way too expensive to
540 ;; check in normal use, and are disabled here.
544 ;; We'll be doing a lot of modular arithmetic.
548 ;;; I'm not sure why I need this FTYPE declaration. Compiled by the
549 ;;; target compiler, it can deduce the return type fine, but without
550 ;;; it, we pay a heavy price in BIGNUM-GCD when compiled by the
551 ;;; cross-compiler. -- CSR, 2004-07-19
554 bignum-factors-of-two))
568 ;;; Multiply a bignum buffer with a fixnum or a digit, storing the
569 ;;; result in another bignum buffer, and without using any
570 ;;; temporaries. Inlined to avoid boxing smallnum if it's actually a
571 ;;; digit. Needed by GCD, should possibly OAOO with
572 ;;; MULTIPLY-BIGNUM-AND-FIXNUM.
575 smallnum res)
590 smallnum
591 carry-digit)
601 res-len)))
603 ;;; Given U and V, return U / V mod 2^32. Implements the algorithm in the
604 ;;; paper, but uses some clever bit-twiddling nicked from Nickle to do it.
620 m))
625 digit-size))
626 do
630 bmod
631 tmp1)))
633 tmp1 tmp1-len
634 u))
648 v))
649 tmp1)))
651 tmp1 tmp1-len
652 u))
654 u-len))
658 ;;; Find D and N such that (LOGAND ALL-ONES-DIGIT (- (* D X) (* N Y))) is 0,
659 ;;; (< 0 N LOWER-ONES-DIGIT) and (< 0 (ABS D) LOWER-ONES-DIGIT).
675 d1 d2
676 n2 n1
677 d2 d1))
687 b))
689 ;;; Allocate a single word bignum that holds fixnum. This is useful when
690 ;;; we are trying to mix fixnum and bignum operands.
695 res))
697 ;; When the larger number is less than this many bignum digits long, revert
698 ;; to old algorithm.
701 ;;; Alternate between k-ary reduction with the help of
702 ;;; REDUCED-RATIO-MOD and digit modulus reduction via DMOD. Once the
703 ;;; arguments get small enough, drop through to BIGNUM-MOD-GCD (since
704 ;;; k-ary reduction can introduce spurious factors, which need to be
705 ;;; filtered out). Reference: Kenneth Weber, "The accelerated integer
706 ;;; GCD algorithm", ACM Transactions on Mathematical Software, volume
707 ;;; 21, number 1, March 1995, epp. 111-122.
711 u0
714 v0
723 n)))
741 buffer-len u-len v-len tmp1-len tmp2-len))
747 factors-of-two)))
751 factors-of-two)))
757 do
766 v v-len
767 tmp1))
771 n tmp1))
774 d tmp2))
781 tmp2 tmp2-len
782 u
784 tmp2-len))))
799 n))
800 factors-of-two)))))
806 ;; While the length difference of A and B is sufficiently large,
807 ;; reduce using MOD (slowish, but it should equalize the sizes of
808 ;; A and B pretty quickly). After that, use the binary GCD
809 ;; algorithm to handle the rest.
820 a))
832 b-buffer len-b))
834 a-buffer
836 factors-of-two)))
838 b-buffer
840 factors-of-two))))
876 res))))
881 a len-a
882 res))))))))
902 increment)))))))
904 \f
905 ;;;; negation
909 ;;; This negates bignum-len digits of bignum, storing the resulting digits into
910 ;;; result (possibly EQ to bignum) and returning whatever end-carry there is.
912 bignum-len
919 (,carry
943 ;;; Fully-normalize is an internal optional. It cause this to always return
944 ;;; a bignum, without any extraneous digits, and it never returns a fixnum.
958 ;;; This assumes bignum is positive; that is, the result of negating it will
959 ;;; stay in the provided allocated bignum.
962 bignum)
971 \f
972 ;;;; shifting
976 ;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
977 ;;; BIGNUM-LDB-BIGNUM-RES. They supply a termination form that references
978 ;;; locals established by this form. Source is the source bignum. Start-digit
979 ;;; is the first digit in source from which we pull bits. Start-pos is the
980 ;;; first bit we want. Res-len-form is the form that computes the length of
981 ;;; the resulting bignum. Termination is a DO termination form with a test and
982 ;;; body. When result is supplied, it is the variable to which this binds a
983 ;;; newly allocated bignum.
984 ;;;
985 ;;; Given start-pos, 1-31 inclusively, of shift, we form the j'th resulting
986 ;;; digit from high bits of the i'th source digit and the start-pos number of
987 ;;; bits from the i+1'th source digit.
989 start-digit
990 start-pos
991 res-len-form
992 termination
1001 ,termination
1007 high-bits-in-first-digit))))))
1011 ;;; First compute the number of whole digits to shift, shifting them by
1012 ;;; skipping them when we start to pick up bits, and the number of bits to
1013 ;;; shift the remaining digits into place. If the number of digits is greater
1014 ;;; than the length of the bignum, then the result is either 0 or -1. If we
1015 ;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
1016 ;;; digits. The last branch handles the general case which uses a macro that a
1017 ;;; couple other routines use. The fifth argument to the macro references
1018 ;;; locals established by the macro.
1038 res)))))
1041 ;; Since a FIXNUM should be big enough to address anything in
1042 ;; memory, including arrays of bits, and since arrays of bits
1043 ;; take up about the same space as corresponding fixnums, there
1044 ;; should be no way that we fall through to this case: any shift
1045 ;; right by a bignum should give zero. But let's check anyway:
1058 ;;; GCD uses this for an in-place shifting operation. This is different enough
1059 ;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
1060 ;;; macro, but they share the basic algorithm. This routine foregoes a first
1061 ;;; test for digits being greater than or equal to bignum-len since that will
1062 ;;; never happen for its uses in GCD. We did fold the last branch into a macro
1063 ;;; since it was duplicated a few times, and the fifth argument to it
1064 ;;; references locals established by the macro.
1082 ;;; This handles shifting a bignum buffer to provide fresh bignum data for some
1083 ;;; internal routines. We know bignum is safe when called with bignum-len.
1084 ;;; First we compute the number of whole digits to shift, shifting them
1085 ;;; starting to store farther along the result bignum. If we shift on a digit
1086 ;;; boundary (that is, n-bits is zero), then we just copy digits. The last
1087 ;;; branch handles the general case.
1101 ;; Left shift by a number too big to be represented as a fixnum
1102 ;; would exceed our memory capacity, since a fixnum is big enough
1103 ;; to index any array, including a bit array.
1113 res))
1115 ;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
1116 ;;; When res comes in non-nil, then this foregoes allocating a result, and it
1117 ;;; normalizes the buffer instead of the would-be allocated result.
1118 ;;;
1119 ;;; We start storing into one digit higher than digits, storing a whole result
1120 ;;; digit from parts of two contiguous digits from bignum. When the loop
1121 ;;; finishes, we store the remaining bits from bignum's first digit in the
1122 ;;; first non-zero result digit, digits. We also grab some left over high
1123 ;;; bits from the last digit of bignum.
1145 remaining-bits)
1147 \f
1148 ;;;; relational operators
1150 ;;; Return T iff bignum is positive.
1155 ;;; This compares two bignums returning -1, 0, or 1, depending on
1156 ;;; whether a is less than, equal to, or greater than b.
1169 (())
1182 \f
1183 ;;;; float conversion
1185 ;;; Make a single or double float with the specified significand,
1186 ;;; exponent and sign.
1197 res
1210 hi
1212 lo)))
1219 exp
1224 ;;; Convert Bignum to a float in the specified Format, rounding to the best
1225 ;;; approximation.
1250 bits
1253 plusp))
1256 bits
1259 plusp))
1260 #!+long-float
1263 bits
1266 plusp))))
1271 ;; Why a SIMPLE-TYPE-ERROR? Well, this is mainly
1272 ;; called by COERCE, which requires an error of
1273 ;; TYPE-ERROR if the conversion can't happen
1274 ;; (except in certain circumstances when we are
1275 ;; coercing to a FUNCTION) -- CSR, 2002-09-18
1279 :expected-type format
1281 exp)))
1284 ;; Round down if round bit is 0.
1287 ;; If only round bit is set, then round to even.
1290 t)
1296 ;; Otherwise, round up.
1299 \f
1300 ;;;; integer length and logbitp/logcount
1334 result
1339 \f
1340 ;;;; logical operations
1342 ;;;; NOT
1353 ;;;; AND
1373 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1374 ;;; is AND, we don't care about any bits longer than a's since its infinite 0
1375 ;;; sign bits will mask the other bits out of b. The result is len-a big.
1385 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1386 ;;; is AND, we just copy any bits longer than a's since its infinite 1 sign
1387 ;;; bits will include any bits from b. The result is len-b big.
1401 ;;;; IOR
1421 ;;; This takes a shorter bignum, a and len-a, that is positive. Because this
1422 ;;; is IOR, we don't care about any bits longer than a's since its infinite
1423 ;;; 0 sign bits will mask the other bits out of b out to len-b. The result
1424 ;;; is len-b long.
1438 ;;; This takes a shorter bignum, a and len-a, that is negative. Because this
1439 ;;; is IOR, we just copy any bits longer than a's since its infinite 1 sign
1440 ;;; bits will include any bits from b. The result is len-b long.
1455 ;;;; XOR
1466 ;;; This takes the shorter of two bignums in a and len-a. Res is len-b
1467 ;;; long. Do the XOR.
1481 \f
1482 ;;;; There used to be a bunch of code to implement "efficient" versions of LDB
1483 ;;;; and DPB here. But it apparently was never used, so it's been deleted.
1484 ;;;; --njf, 2007-02-04
1485 \f
1486 ;;;; TRUNCATE
1488 ;;; This is the original sketch of the algorithm from which I implemented this
1489 ;;; TRUNCATE, assuming both operands are bignums. I should modify this to work
1490 ;;; with the documentation on my functions, as a general introduction. I've
1491 ;;; left this here just in case someone needs it in the future. Don't look at
1492 ;;; this unless reading the functions' comments leaves you at a loss. Remember
1493 ;;; this comes from Knuth, so the book might give you the right general
1494 ;;; overview.
1495 ;;;
1496 ;;; (truncate x y):
1497 ;;;
1498 ;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
1499 ;;;
1500 ;;; Make x and y positive, copying x if it is already positive.
1501 ;;;
1502 ;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
1503 ;;; digit)
1504 ;;; Just do most sig digit to determine how much to shift whole number.
1505 ;;; Shift x this much too.
1506 ;;; Remember this initial shift count.
1507 ;;;
1508 ;;; Allocate q to be len-x minus len-y quantity plus 1.
1509 ;;;
1510 ;;; i = last digit of x.
1511 ;;; k = last digit of q.
1512 ;;;
1513 ;;; LOOP
1514 ;;;
1515 ;;; j = last digit of y.
1516 ;;;
1517 ;;; compute guess.
1518 ;;; if x[i] = y[j] then g = (1- (ash 1 digit-size))
1519 ;;; else g = x[i]x[i-1]/y[j].
1520 ;;;
1521 ;;; check guess.
1522 ;;; %UNSIGNED-MULTIPLY returns b and c defined below.
1523 ;;; a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
1524 ;;; Use %UNSIGNED-MULTIPLY taking low-order result.
1525 ;;; b = (logand (ash (* g y[j-1]) (- digit-size)) (1- (ash 1 digit-size))).
1526 ;;; c = (logand (* g y[j-1]) (1- (ash 1 digit-size))).
1527 ;;; if a < b, okay.
1528 ;;; if a > b, guess is too high
1529 ;;; g = g - 1; go back to "check guess".
1530 ;;; if a = b and c > x[i-2], guess is too high
1531 ;;; g = g - 1; go back to "check guess".
1532 ;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
1533 ;;; SAME FOR A, B, AND C.
1534 ;;;
1535 ;;; Subtract g * y from x[i - len-y+1]..x[i]. See paper for doing this in step.
1536 ;;; If x[i] < 0, guess is screwed up.
1537 ;;; negative g, then add 1
1538 ;;; zero or positive g, then subtract 1
1539 ;;; AND add y back into x[len-y+1..i].
1540 ;;;
1541 ;;; q[k] = g.
1542 ;;; i = i - 1.
1543 ;;; k = k - 1.
1544 ;;;
1545 ;;; If k>=0, goto LOOP.
1546 ;;;
1547 ;;; Now quotient is good, but remainder is not.
1548 ;;; Shift x right by saved initial left shifting count.
1549 ;;;
1550 ;;; Check quotient and remainder signs.
1551 ;;; x pos y pos --> q pos r pos
1552 ;;; x pos y neg --> q neg r pos
1553 ;;; x neg y pos --> q neg r neg
1554 ;;; x neg y neg --> q pos r neg
1555 ;;;
1556 ;;; Normalize quotient and remainder. Cons result if necessary.
1559 ;;; This used to be split into multiple functions, which shared state
1560 ;;; in special variables *TRUNCATE-X* and *TRUNCATE-Y*. Having so many
1561 ;;; special variable accesses in tight inner loops was having a large
1562 ;;; effect on performance, so the helper functions have now been
1563 ;;; refactored into local functions and the special variables into
1564 ;;; lexicals. There was also a lot of boxing and unboxing of
1565 ;;; (UNSIGNED-BYTE 32)'s going on, which this refactoring
1566 ;;; eliminated. This improves the performance on some CL-BENCH tests
1567 ;;; by up to 50%, which is probably signigicant enough to justify the
1568 ;;; reduction in readability that was introduced. --JES, 2004-08-07
1573 ;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
1574 ;;; fixes up the quotient and remainder with respect to sign and
1575 ;;; normalization.
1576 ;;;
1577 ;;; We don't have to worry about shifting Y to make its most
1578 ;;; significant digit sufficiently large for %FLOOR to return
1579 ;;; digit-size quantities for the q-digit and r-digit. If Y is
1580 ;;; a single digit bignum, it is already large enough for
1581 ;;; %FLOOR. That is, it has some bits on pretty high in the
1582 ;;; digit.
1588 ;; Y is a power of two.
1589 ;; SHIFT-RIGHT-UNALIGNED won't do the right thing
1590 ;; with a shift count of 0 or -1, so special case this.
1594 ;; We could probably get away with (VALUES X 0)
1595 ;; here, but it's not clear that some of the
1596 ;; normalization logic further down would avoid
1597 ;; mutilating X. Just go ahead and cons, consing's
1598 ;; cheap.
1607 res)
1608 res)
1623 ;;; This returns a guess for the next division step. Y1 is the
1624 ;;; highest y digit, and y2 is the second to highest y
1625 ;;; digit. The x... variables are the three highest x digits
1626 ;;; for the next division step.
1627 ;;;
1628 ;;; From Knuth, our guess is either all ones or x-i and x-i-1
1629 ;;; divided by y1, depending on whether x-i and y1 are the
1630 ;;; same. We test this guess by determining whether guess*y2
1631 ;;; is greater than the three high digits of x minus guess*y1
1632 ;;; shifted left one digit:
1633 ;;; ------------------------------
1634 ;;; | x-i | x-i-1 | x-i-2 |
1635 ;;; ------------------------------
1636 ;;; ------------------------------
1637 ;;; - | g*y1 high | g*y1 low | 0 |
1638 ;;; ------------------------------
1639 ;;; ... < guess*y2 ???
1640 ;;; If guess*y2 is greater, then we decrement our guess by one
1641 ;;; and try again. This returns a guess that is either
1642 ;;; correct or one too large.
1646 all-ones-digit
1662 ;; Supplying borrow of 1 means there was no
1663 ;; borrow, and we know x-i-2 minus 0 requires
1664 ;; no borrow.
1666 high-guess*y1
1667 borrow)))
1671 middle-digit)
1675 x-i-2))))
1678 ;;; Divide TRUNCATE-X by TRUNCATE-Y, returning the quotient
1679 ;;; and destructively modifying TRUNCATE-X so that it holds
1680 ;;; the remainder.
1681 ;;;
1682 ;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
1683 ;;;
1684 ;;; TRUNCATE-X definitely has at least three digits, and it has one
1685 ;;; more than TRUNCATE-Y. This keeps i, i-1, i-2, and low-x-digit
1686 ;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
1690 ;; Add one for extra sign digit in case high bit is on.
1704 ;; This modifies TRUNCATE-X. Must access
1705 ;; elements each pass.
1710 len-y low-x-digit))
1715 q))
1716 ;;; This takes a digit guess, multiplies it by TRUNCATE-Y for a
1717 ;;; result one greater in length than LEN-Y, and subtracts this result
1718 ;;; from TRUNCATE-X. LOW-X-DIGIT is the first digit of X to start
1719 ;;; the subtraction, and we know X is long enough to subtract a LEN-Y
1720 ;;; plus one length bignum from it. Next we check the result of the
1721 ;;; subtraction, and if the high digit in X became negative, then our
1722 ;;; guess was one too big. In this case, return one less than GUESS
1723 ;;; passed in, and add one value of Y back into X to account for
1724 ;;; subtracting one too many. Knuth shows that the guess is wrong on
1725 ;;; the order of 3/b, where b is the base (2 to the digit-size power)
1726 ;;; -- pretty rarely.
1736 ;; Multiply guess and divisor, subtracting from dividend
1737 ;; simultaneously.
1740 (%multiply-and-add guess
1742 carry-digit)
1747 low-digit
1748 borrow)
1756 carry-digit borrow))
1757 ;; See whether guess is off by one, adding one
1758 ;; Y back in if necessary.
1760 guess)
1762 ;; If subtraction has negative result, add one
1763 ;; divisor value back in. The guess was one too
1764 ;; large in magnitude.
1771 carry)
1780 ;;; This returns the amount to shift y to place a one in the
1781 ;;; second highest bit. Y must be positive. If the last digit
1782 ;;; of y is zero, then y has a one in the previous digit's
1783 ;;; sign bit, so we know it will take one less than digit-size
1784 ;;; to get a one where we want. Otherwise, we count how many
1785 ;;; right shifts it takes to get zero; subtracting this value
1786 ;;; from digit-size tells us how many high zeros there are
1787 ;;; which is one more than the shift amount sought.
1788 ;;;
1789 ;;; Note: This is exactly the same as one less than the
1790 ;;; integer-length of the last digit subtracted from the
1791 ;;; digit-size.
1792 ;;;
1793 ;;; We shift y to make it sufficiently large that doing the
1794 ;;; 2*digit-size by digit-size %FLOOR calls ensures the quotient and
1795 ;;; remainder fit in digit-size.
1802 ;;; Stores two bignums into the truncation bignum buffers,
1803 ;;; shifting them on the way in. This assumes x and y are
1804 ;;; positive and at least two in length, and it assumes
1805 ;;; truncate-x and truncate-y are one digit longer than x and
1806 ;;; y.
1815 truncate-x)
1818 ;;; Divide X by Y returning the quotient and remainder. In the
1819 ;;; general case, we shift Y to set up for the algorithm, and we
1820 ;;; use two buffers to save consing intermediate values. X gets
1821 ;;; destructively modified to become the remainder, and we have
1822 ;;; to shift it to account for the initial Y shift. After we
1823 ;;; multiple bind q and r, we first fix up the signs and then
1824 ;;; return the normalized results.
1845 len-y y-shift)
1847 len-y)
1848 ;; Now that RETURN-QUOTIENT-LEAVING-REMAINDER
1849 ;; has executed, we just tidy up the remainder
1850 ;; (in TRUNCATE-X) and return it.
1859 truncate-x 0 y-shift len-y
1863 y-shift))
1865 res))))))))
1873 quotient
1876 rem
1879 \f
1880 ;;;; There used to be a pile of code for implementing division for bignum digits
1881 ;;;; for machines that don't have a 2*digit-size by digit-size divide instruction.
1882 ;;;; This happens to be most machines, but all the SBCL ports seem to be content
1883 ;;;; to implement SB-BIGNUM:%FLOOR as a VOP rather than using the code here.
1884 ;;;; So it's been deleted. --njf, 2007-02-04
1885 \f
1886 ;;;; general utilities
1888 ;;; Internal in-place operations use this to fixup remaining digits in the
1889 ;;; incoming data, such as in-place shifting. This is basically the same as
1890 ;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
1891 ;;; instead of shrinking the bignum.
1905 len)
1907 ;;; This drops the last digit if it is unnecessary sign information. It repeats
1908 ;;; this as needed, possibly ending with a fixnum. If the resulting length from
1909 ;;; shrinking is one, see whether our one word is a fixnum. Shift the possible
1910 ;;; fixnum bits completely out of the word, and compare this with shifting the
1911 ;;; sign bit all the way through. If the bits are all 1's or 0's in both words,
1912 ;;; then there are just sign bits between the fixnum bits and the sign bit. If
1913 ;;; we do have a fixnum, shift it over for the two low-tag bits.
1927 result))
1928 result)))
1930 ;;; This drops the last digit if it is unnecessary sign information. It
1931 ;;; repeats this as needed, possibly ending with a fixnum magnitude but never
1932 ;;; returning a fixnum.
1941 result))
1942 \f
1943 ;;;; hashing
1945 ;;; the bignum case of the SXHASH function
1956 result))