home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / pcl / src-16f.lha / code / bignum.lisp < prev    next >
Encoding:
Text File  |  1992-05-30  |  89.3 KB  |  2,564 lines

  1. ;;; -*- Mode: completion; Log: code.log; Package: bignum -*-
  2. ;;;
  3. ;;; **********************************************************************
  4. ;;; This code was written as part of the CMU Common Lisp project at
  5. ;;; Carnegie Mellon University, and has been placed in the public domain.
  6. ;;; If you want to use this code or any part of CMU Common Lisp, please contact
  7. ;;; Scott Fahlman or slisp-group@cs.cmu.edu.
  8. ;;;
  9. (ext:file-comment
  10.   "$Header: bignum.lisp,v 1.19 91/06/12 17:24:18 chiles Exp $")
  11. ;;;
  12. ;;; **********************************************************************
  13. ;;;
  14. ;;; This file contains code to implement bignum support.
  15. ;;;
  16.  
  17. (in-package "BIGNUM")
  18. (use-package "KERNEL")
  19.  
  20. ;;; These symbols define the interface to the number code.
  21.  
  22. (export '(add-bignums multiply-bignums negate-bignum subtract-bignum
  23.       multiply-bignum-and-fixnum multiply-fixnums
  24.       bignum-ashift-right bignum-ashift-left bignum-gcd
  25.       bignum-to-float float-bignum-ratio bignum-integer-length
  26.       bignum-logical-and bignum-logical-ior bignum-logical-xor
  27.       bignum-logical-not bignum-load-byte bignum-deposit-byte
  28.       bignum-truncate bignum-plus-p bignum-compare make-small-bignum
  29.       bignum-logcount))
  30.  
  31. ;;; These symbols define the interface to the compiler.
  32.  
  33. (export '(bignum-type bignum-element-type bignum-index %allocate-bignum
  34.       %bignum-length %bignum-set-length %bignum-ref %bignum-set
  35.       %digit-0-or-plusp %add-with-carry %subtract-with-borrow
  36.       %multiply-and-add %multiply %lognot %logand %logior %logxor
  37.       %fixnum-to-digit %floor %fixnum-digit-with-correct-sign %ashl
  38.       %ashr %digit-logical-shift-right))
  39.  
  40.  
  41.  
  42. ;;;; Notes.
  43.  
  44. ;;; The following interfaces will either be assembler routines or code sequences
  45. ;;; expanded into the code as basic bignum operations:
  46. ;;;    General:
  47. ;;;       %BIGNUM-LENGTH
  48. ;;;       %ALLOCATE-BIGNUM
  49. ;;;       %BIGNUM-REF
  50. ;;;       %NORMALIZE-BIGNUM
  51. ;;;       %BIGNUM-SET-LENGTH
  52. ;;;       %FIXNUM-DIGIT-WITH-CORRECT-SIGN
  53. ;;;       %SIGN-DIGIT
  54. ;;;      %ASHR
  55. ;;;       %ASHL
  56. ;;;       %BIGNUM-0-OR-PLUSP
  57. ;;;       %DIGIT-LOGICAL-SHIFT-RIGHT
  58. ;;;    General (May not exist when done due to sole use in %-routines.)
  59. ;;;       %DIGIT-0-OR-PLUSP
  60. ;;;    Addition:
  61. ;;;       %ADD-WITH-CARRY
  62. ;;;    Subtraction:
  63. ;;;       %SUBTRACT-WITH-BORROW
  64. ;;;    Multiplication
  65. ;;;       %MULTIPLY
  66. ;;;    Negation
  67. ;;;       %LOGNOT
  68. ;;;    Shifting (in place)
  69. ;;;       %NORMALIZE-BIGNUM-BUFFER
  70. ;;;    GCD/Relational operators:
  71. ;;;       %DIGIT-COMPARE
  72. ;;;       %DIGIT-GREATER
  73. ;;;    Relational operators:
  74. ;;;       %LOGAND
  75. ;;;       %LOGIOR
  76. ;;;       %LOGXOR
  77. ;;;    LDB
  78. ;;;       %FIXNUM-TO-DIGIT
  79. ;;;    TRUNCATE
  80. ;;;       %FLOOR
  81. ;;;
  82. ;;;
  83. ;;; Note: The floating routines know about the float representation.
  84. ;;;
  85. ;;; PROBLEM 1:
  86. ;;; There might be a problem with various LET's and parameters that take a
  87. ;;; digit value.  We need to write these so those things stay in 32-bit
  88. ;;; registers and number stack slots.  I bind locals to these values, and I
  89. ;;; use function on them -- ZEROP, ASH, etc.
  90. ;;;
  91. ;;; PROBLEM 2:
  92. ;;; In shifting and byte operations, I use masks and logical operations that
  93. ;;; could result in intermediate bignums.  This is hidden by the current system,
  94. ;;; but I may need to write these in a way that keeps these masks and logical
  95. ;;; operations from diving into the Lisp level bignum code.
  96. ;;;
  97. ;;; To do:
  98. ;;;    fixnums
  99. ;;;       logior, logxor, logand
  100. ;;;       depending on relationals, < (twice) and <= (twice)
  101. ;;;          or write compare thing (twice).
  102. ;;;       LDB on fixnum with bignum result.
  103. ;;;       DPB on fixnum with bignum result.
  104. ;;;       TRUNCATE returns zero or one as one value and fixnum or minus fixnum
  105. ;;;          for the other value when given (truncate fixnum bignum).
  106. ;;;          Returns (truncate bignum fixnum) otherwise.
  107. ;;;       addition
  108. ;;;       subtraction (twice)
  109. ;;;       multiply
  110. ;;;       GCD
  111. ;;;    write MASK-FIELD and DEPOSIT-FIELD in terms of logical operations.
  112. ;;;    DIVIDE
  113. ;;;       IF (/ x y) with bignums:
  114. ;;;          do the truncate, and if rem is 0, return quotient.
  115. ;;;          if rem is non-0
  116. ;;;         gcd of x and y.
  117. ;;;         "truncate" each by gcd, ignoring remainder 0.
  118. ;;;         form ratio of each result, bottom is positive.
  119. ;;;
  120.  
  121.  
  122.  
  123. ;;;; What's a bignum?
  124.  
  125. (eval-when (compile load eval) ;Necessary for DEFTYPE.
  126.  
  127. (defconstant digit-size 32)
  128.  
  129. (defconstant maximum-bignum-length (1- (ash 1 (- 32 vm:type-bits))))
  130.  
  131. ) ;eval-when
  132.  
  133.  
  134.  
  135. ;;;; Internal inline routines.
  136.  
  137. ;;; %ALLOCATE-BIGNUM must zero all elements.
  138. ;;;
  139. (defun %allocate-bignum (length)
  140.   (declare (type bignum-index length))
  141.   (%allocate-bignum length))
  142.  
  143. ;;; Extract the length of the bignum.
  144. ;;; 
  145. (defun %bignum-length (bignum)
  146.   (declare (type bignum-type bignum))
  147.   (%bignum-length bignum))
  148.  
  149. ;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
  150. ;;; to be able to return 32 bits somewhere no one looks for real objects.
  151. ;;;
  152. (defun %bignum-ref (bignum i)
  153.   (declare (type bignum-type bignum)
  154.        (type bignum-index i))
  155.   (%bignum-ref bignum i))
  156. ;;;
  157. (defun %bignum-set (bignum i value)
  158.   (declare (type bignum-type bignum)
  159.        (type bignum-index i)
  160.        (type bignum-element-type value))
  161.   (%bignum-set bignum i value))
  162. ;;;
  163. (defsetf %bignum-ref %bignum-set)
  164.  
  165. ;;; Return T if digit is positive, or NIL if negative.
  166. ;;;
  167. (defun %digit-0-or-plusp (digit)
  168.   (declare (type bignum-element-type digit))
  169.   (not (logbitp (1- digit-size) digit)))
  170.  
  171. (proclaim '(inline %bignum-0-or-plusp))
  172. (defun %bignum-0-or-plusp (bignum len)
  173.   (declare (type bignum-type bignum)
  174.        (type bignum-index len))
  175.   (%digit-0-or-plusp (%bignum-ref bignum (1- len))))
  176.  
  177. ;;; %ADD-WITH-CARRY -- Internal.
  178. ;;;
  179. ;;; This should be in assembler, and should not cons intermediate results.  It
  180. ;;; returns a 32bit digit and a carry resulting from adding together a, b, and
  181. ;;; an incoming carry.
  182. ;;;
  183. (defun %add-with-carry (a b carry)
  184.   (declare (type bignum-element-type a b)
  185.        (type (mod 2) carry))
  186.   (%add-with-carry a b carry))
  187.  
  188. ;;; %SUBTRACT-WITH-BORROW -- Internal.
  189. ;;;
  190. ;;; This should be in assembler, and should not cons intermediate results.  It
  191. ;;; returns a 32bit digit and a borrow resulting from subtracting b from a, and
  192. ;;; subtracting a possible incoming borrow.
  193. ;;;
  194. ;;; We really do:  a - b - 1 + borrow, where borrow is either 0 or 1.
  195. ;;; 
  196. (defun %subtract-with-borrow (a b borrow)
  197.   (declare (type bignum-element-type a b)
  198.        (type (mod 2) borrow))
  199.   (%subtract-with-borrow a b borrow))
  200.  
  201. ;;; %MULTIPLY -- Internal.
  202. ;;;
  203. ;;; This multiplies two digit-size (32-bit) numbers, returning a 64-bit result
  204. ;;; split into two 32-bit quantities.
  205. ;;;
  206. (defun %multiply (x y)
  207.   (declare (type bignum-element-type x y))
  208.   (%multiply x y))
  209.  
  210. ;;; %MULTIPLY-AND-ADD  --  Internal.
  211. ;;;
  212. ;;; This multiplies x-digit and y-digit, producing high and low digits
  213. ;;; manifesting the result.  Then it adds the low digit, res-digit, and
  214. ;;; carry-in-digit.  Any carries (note, you still have to add two digits at a
  215. ;;; time possibly producing two carries) from adding these three digits get
  216. ;;; added to the high digit from the multiply, producing the next carry digit.
  217. ;;; Res-digit is optional since two uses of this primitive multiplies a single
  218. ;;; digit bignum by a multiple digit bignum, and in this situation there is no
  219. ;;; need for a result buffer accumulating partial results which is where the
  220. ;;; res-digit comes from.
  221. ;;;
  222. (defun %multiply-and-add (x-digit y-digit carry-in-digit &optional (res-digit 0))
  223.   (declare (type bignum-element-type x-digit y-digit res-digit carry-in-digit))
  224.   (%multiply-and-add x-digit y-digit carry-in-digit res-digit))
  225.  
  226. ;;; %LOGNOT -- Internal.
  227. ;;;
  228. (defun %lognot (digit)
  229.   (declare (type bignum-element-type digit))
  230.   (%lognot digit))
  231.  
  232. ;;; %LOGAND -- Internal.
  233. ;;; %LOGIOR -- Internal.
  234. ;;; %LOGXOR -- Internal.
  235. ;;;
  236. ;;; Do the 32bit unsigned op.
  237. ;;;
  238. (proclaim '(inline %logand %logior %logxor))
  239. (defun %logand (a b)
  240.   (declare (type bignum-element-type a b))
  241.   (logand a b))
  242. (defun %logior (a b)
  243.   (declare (type bignum-element-type a b))
  244.   (logior a b))
  245. (defun %logxor (a b)
  246.   (declare (type bignum-element-type a b))
  247.   (logxor a b))
  248.  
  249. ;;; %FIXNUM-TO-DIGIT -- Internal.
  250. ;;;
  251. ;;; This takes a fixnum and sets it up as an unsigned 32-bit quantity.  In
  252. ;;; the new system this will mean shifting it right two bits.
  253. ;;;
  254. (defun %fixnum-to-digit (x)
  255.   (declare (fixnum x))
  256.   (logand x (1- (ash 1 digit-size))))
  257.  
  258. #-32x16-divide
  259. ;;; %FLOOR -- Internal.
  260. ;;;
  261. ;;; This takes three digits and returns the FLOOR'ed result of dividing the
  262. ;;; first two as a 64-bit integer by the third.
  263. ;;;
  264. ;;; DO WEIRD let AND setq STUFF TO SLIME THE COMPILER INTO ALLOWING THE %FLOOR
  265. ;;; TRANSFORM TO EXPAND INTO PSEUDO-ASSEMBLER FOR WHICH THE COMPILER CAN LATER
  266. ;;; CORRECTLY ALLOCATE REGISTERS.
  267. ;;;
  268. (defun %floor (a b c)
  269.   (let ((a a) (b b) (c c))
  270.     (declare (type bignum-element-type a b c))
  271.     (setq a a b b c c)
  272.     (%floor a b c)))
  273.  
  274.  
  275. ;;; %FIXNUM-DIGIT-WITH-CORRECT-SIGN -- Internal.
  276. ;;;
  277. ;;; Convert the digit to a regular integer assuming that the digit is signed.
  278. ;;;
  279. (defun %fixnum-digit-with-correct-sign (digit)
  280.   (declare (type bignum-element-type digit))
  281.   (if (logbitp (1- digit-size) digit)
  282.       (logior digit (ash -1 digit-size))
  283.       digit))
  284.  
  285. ;;; %ASHR -- Internal.
  286. ;;;
  287. ;;; Do an arithmetic shift right of data even though bignum-element-type is
  288. ;;; unsigned.
  289. ;;;
  290. (defun %ashr (data count)
  291.   (declare (type bignum-element-type data)
  292.        (type (mod 32) count))
  293.   (%ashr data count))
  294.  
  295. ;;; %ASHL -- Internal.
  296. ;;;
  297. ;;; This takes a 32-bit quantity and shifts it to the left, returning a 32-bit
  298. ;;; quantity.
  299. (defun %ashl (data count)
  300.   (declare (type bignum-element-type data)
  301.        (type (mod 32) count))
  302.   (%ashl data count))
  303.  
  304. ;;; %DIGIT-LOGICAL-SHIFT-RIGHT  --  Internal
  305. ;;;
  306. ;;;    Do an unsigned (logical) right shift of a digit by Count.
  307. ;;;
  308. (defun %digit-logical-shift-right (data count)
  309.   (declare (type bignum-element-type data)
  310.        (type (mod 32) count))
  311.   (%digit-logical-shift-right data count))
  312.  
  313.  
  314. ;;; %BIGNUM-SET-LENGTH -- Internal.
  315. ;;;
  316. ;;; Change the length of bignum to be newlen.  Newlen must be the same or
  317. ;;; smaller than the old length, and any elements beyond newlen must be zeroed.
  318. ;;; 
  319. (defun %bignum-set-length (bignum newlen)
  320.   (declare (type bignum-type bignum)
  321.        (type bignum-index newlen))
  322.   (%bignum-set-length bignum newlen))
  323.  
  324. ;;; %SIGN-DIGIT -- Internal.
  325. ;;;
  326. ;;; This returns 0 or "-1" depending on whether the bignum is positive.  This
  327. ;;; is suitable for infinite sign extension to complete additions,
  328. ;;; subtractions, negations, etc.  This cannot return a -1 represented as
  329. ;;; a negative fixnum since it would then have to low zeros.
  330. ;;;
  331. (proclaim '(inline %sign-digit))
  332. (defun %sign-digit (bignum len)
  333.   (declare (type bignum-type bignum)
  334.        (type bignum-index len))
  335.   (%ashr (%bignum-ref bignum (1- len)) (1- digit-size)))
  336.  
  337.  
  338.  
  339. ;;; %DIGIT-COMPARE and %DIGIT-GREATER -- Internal.
  340. ;;;
  341. ;;; These take two 32 bit quantities and compare or contrast them without
  342. ;;; wasting time with incorrect type checking.
  343. ;;;
  344. (proclaim '(inline %digit-compare %digit-greater))
  345. (defun %digit-compare (x y)
  346.   (= x y))
  347. ;;;
  348. (defun %digit-greater (x y)
  349.   (> x y))
  350.  
  351.  
  352. (proclaim '(optimize (speed 3) (safety 0)))
  353.  
  354.  
  355. ;;;; Addition.
  356.  
  357. (defun add-bignums (a b)
  358.   (declare (type bignum-type a b))
  359.   (let ((len-a (%bignum-length a))
  360.     (len-b (%bignum-length b)))
  361.     (declare (type bignum-index len-a len-b))
  362.     (multiple-value-bind (a len-a b len-b)
  363.              (if (> len-a len-b)
  364.                  (values a len-a b len-b)
  365.                  (values b len-b a len-a))
  366.       (declare (type bignum-type a b)
  367.            (type bignum-index len-a len-b))
  368.       (let* ((len-res (1+ len-a))
  369.          (res (%allocate-bignum len-res))
  370.          (carry 0))
  371.     (declare (type bignum-index len-res)
  372.          (type bignum-type res)
  373.          (type (mod 2) carry))
  374.     (dotimes (i len-b)
  375.       (declare (type bignum-index i))
  376.       (multiple-value-bind
  377.           (v k)
  378.           (%add-with-carry (%bignum-ref a i) (%bignum-ref b i) carry)
  379.         (declare (type bignum-element-type v)
  380.              (type (mod 2) k))
  381.         (setf (%bignum-ref res i) v)
  382.         (setf carry k)))
  383.     (if (/= len-a len-b)
  384.         (finish-add a res carry (%sign-digit b len-b) len-b len-a)
  385.         (setf (%bignum-ref res len-a)
  386.           (%add-with-carry (%sign-digit a len-a)
  387.                    (%sign-digit b len-b)
  388.                    carry)))
  389.     (%normalize-bignum res len-res)))))
  390.  
  391. ;;; FINISH-ADD -- Internal.
  392. ;;;
  393. ;;; This takes the longer of two bignums and propagates the carry through its
  394. ;;; remaining high order digits.
  395. ;;;
  396. (defun finish-add (a res carry sign-digit-b start end)
  397.   (declare (type bignum-type a res)
  398.        (type (mod 2) carry)
  399.        (type bignum-element-type sign-digit-b)
  400.        (type bignum-index start end))
  401.   (do ((i start (1+ i)))
  402.       ((= i end)
  403.        (setf (%bignum-ref res end)
  404.          (%add-with-carry (%sign-digit a end) sign-digit-b carry)))
  405.     (declare (type bignum-index i))
  406.     (multiple-value-bind (v k)
  407.              (%add-with-carry (%bignum-ref a i) sign-digit-b carry)
  408.       (setf (%bignum-ref res i) v)
  409.       (setf carry k)))
  410.   (ext:undefined-value))
  411.  
  412.  
  413. ;;;; Subtraction.
  414.  
  415. (eval-when (compile eval)
  416.  
  417. ;;; SUBTRACT-BIGNUM-LOOP -- Internal.
  418. ;;;
  419. ;;; This subtracts b from a plugging result into res.  Return-fun is the
  420. ;;; function to call that fixes up the result returning any useful values, such
  421. ;;; as the result.  This macro may evaluate its arguments more than once.
  422. ;;;
  423. (defmacro subtract-bignum-loop (a len-a b len-b res len-res return-fun)
  424.   (let ((borrow (gensym))
  425.     (a-digit (gensym))
  426.     (a-sign (gensym))
  427.     (b-digit (gensym))
  428.     (b-sign (gensym))
  429.     (i (gensym))
  430.     (v (gensym))
  431.     (k (gensym)))
  432.     `(let* ((,borrow 1)
  433.         (,a-sign (%sign-digit ,a ,len-a))
  434.         (,b-sign (%sign-digit ,b ,len-b)))
  435.        (declare (type bignum-element-type ,a-sign ,b-sign))
  436.        (dotimes (,i ,len-res)
  437.      (declare (type bignum-index ,i))
  438.      (let ((,a-digit (if (< ,i ,len-a) (%bignum-ref ,a ,i) ,a-sign))
  439.            (,b-digit (if (< ,i ,len-b) (%bignum-ref ,b ,i) ,b-sign)))
  440.        (declare (type bignum-element-type ,a-digit ,b-digit))
  441.        (multiple-value-bind
  442.            (,v ,k)
  443.            (%subtract-with-borrow ,a-digit ,b-digit ,borrow)
  444.          (setf (%bignum-ref ,res ,i) ,v)
  445.          (setf ,borrow ,k))))
  446.        (,return-fun ,res ,len-res))))
  447.  
  448. ) ;EVAL-WHEN
  449.  
  450. (defun subtract-bignum (a b)
  451.   (declare (type bignum-type a b))
  452.   (let* ((len-a (%bignum-length a))
  453.      (len-b (%bignum-length b))
  454.      (len-res (1+ (max len-a len-b)))
  455.      (res (%allocate-bignum len-res)))
  456.     (declare (type bignum-index len-a len-b len-res)) ;Test len-res for bounds?
  457.     (subtract-bignum-loop a len-a b len-b res len-res %normalize-bignum)))
  458.  
  459. ;;; SUBTRACT-BIGNUM-BUFFERS -- Internal.
  460. ;;;
  461. ;;; Operations requiring a subtraction without the overhead of intermediate
  462. ;;; results, such as GCD, use this.  It assumes Result is big enough for the
  463. ;;; result.
  464. ;;;
  465. (defun subtract-bignum-buffers (a len-a b len-b result)
  466.   (declare (type bignum-type a b)
  467.        (type bignum-index len-a len-b))
  468.   (let ((len-res (max len-a len-b)))
  469.     (subtract-bignum-loop a len-a b len-b result len-res
  470.               %normalize-bignum-buffer)))
  471.  
  472.  
  473.  
  474. ;;;; Multiplication.
  475.  
  476. (defun multiply-bignums (a b)
  477.   (declare (type bignum-type a b))
  478.   (let* ((a-plusp (%bignum-0-or-plusp a (%bignum-length a)))
  479.      (b-plusp (%bignum-0-or-plusp b (%bignum-length b)))
  480.      (a (if a-plusp a (negate-bignum a)))
  481.      (b (if b-plusp b (negate-bignum b)))
  482.      (len-a (%bignum-length a))
  483.      (len-b (%bignum-length b))
  484.      (len-res (+ len-a len-b))
  485.      (res (%allocate-bignum len-res))
  486.      (negate-res (not (eq a-plusp b-plusp))))
  487.     (declare (type bignum-index len-a len-b len-res))
  488.     (dotimes (i len-a)
  489.       (declare (type bignum-index i))
  490.       (let ((carry-digit 0)
  491.         (x (%bignum-ref a i))
  492.         (k i))
  493.     (declare (type bignum-index k)
  494.          (type bignum-element-type carry-digit x))
  495.     (dotimes (j len-b)
  496.       (multiple-value-bind (big-carry res-digit)
  497.                    (%multiply-and-add x (%bignum-ref b j)
  498.                           (%bignum-ref res k)
  499.                           carry-digit)
  500.         (declare (type bignum-element-type big-carry res-digit))
  501.         (setf (%bignum-ref res k) res-digit)
  502.         (setf carry-digit big-carry)
  503.         (incf k)))
  504.     (setf (%bignum-ref res k) carry-digit)))
  505.     (when negate-res (negate-bignum-in-place res))
  506.     (%normalize-bignum res len-res)))
  507.  
  508. (defun multiply-bignum-and-fixnum (bignum fixnum)
  509.   (declare (type bignum-type bignum) (type fixnum fixnum))
  510.   (let* ((bignum-plus-p (%bignum-0-or-plusp bignum (%bignum-length bignum)))
  511.      (fixnum-plus-p (not (minusp fixnum)))
  512.      (bignum (if bignum-plus-p bignum (negate-bignum bignum)))
  513.      (bignum-len (%bignum-length bignum))
  514.      (fixnum (%fixnum-to-digit (if fixnum-plus-p fixnum (- fixnum))))
  515.      (result (%allocate-bignum (1+ bignum-len)))
  516.      (carry-digit 0))
  517.     (declare (type bignum-type bignum result)
  518.          (type bignum-index bignum-len)
  519.          (type bignum-element-type fixnum carry-digit))
  520.     (dotimes (index bignum-len)
  521.       (declare (type bignum-index index))
  522.       (multiple-value-bind
  523.       (next-digit low)
  524.       (%multiply-and-add (%bignum-ref bignum index) fixnum carry-digit)
  525.     (declare (type bignum-element-type next-digit low))
  526.     (setf carry-digit next-digit)
  527.     (setf (%bignum-ref result index) low)))
  528.     (setf (%bignum-ref result bignum-len) carry-digit)
  529.     (unless (eq bignum-plus-p fixnum-plus-p)
  530.       (negate-bignum-in-place result))
  531.     (%normalize-bignum result (1+ bignum-len))))
  532.  
  533. (defun multiply-fixnums (a b)
  534.   (declare (fixnum a b))
  535.   (let* ((a-minusp (minusp a))
  536.      (b-minusp (minusp b)))
  537.     (multiple-value-bind (high low)
  538.              (%multiply (%fixnum-to-digit (if a-minusp (- a) a))
  539.                     (%fixnum-to-digit (if b-minusp (- b) b)))
  540.       (declare (type bignum-element-type high low))
  541.       (if (and (zerop high)
  542.            (%digit-0-or-plusp low))
  543.       (let ((low (ext:truly-the (unsigned-byte 31)
  544.                     (%fixnum-digit-with-correct-sign low))))
  545.         (if (eq a-minusp b-minusp)
  546.         low
  547.         (- low)))
  548.       (let ((res (%allocate-bignum 2)))
  549.         (%bignum-set res 0 low)
  550.         (%bignum-set res 1 high)
  551.         (unless (eq a-minusp b-minusp) (negate-bignum-in-place res))
  552.         (%normalize-bignum res 2))))))
  553.  
  554.  
  555.  
  556. ;;;; BIGNUM-REPLACE and WITH-BIGNUM-BUFFERS.
  557.  
  558. (eval-when (compile eval)
  559.  
  560. ;;; BIGNUM-REPLACE -- Internal.
  561. ;;;
  562. (defmacro bignum-replace (dest src &key (start1 '0) end1 (start2 '0) end2
  563.                    from-end)
  564.   (ext:once-only ((n-dest dest)
  565.           (n-src src))
  566.     (let ((n-start1 (gensym))
  567.       (n-end1 (gensym))
  568.       (n-start2 (gensym))
  569.       (n-end2 (gensym))
  570.       (i1 (gensym))
  571.       (i2 (gensym))
  572.       (end1 (or end1 `(%bignum-length ,n-dest)))
  573.       (end2 (or end2 `(%bignum-length ,n-src))))
  574.       (if from-end
  575.       `(let ((,n-start1 ,start1)
  576.          (,n-start2 ,start2))
  577.          (do ((,i1 (1- ,end1) (1- ,i1))
  578.           (,i2 (1- ,end2) (1- ,i2)))
  579.          ((or (< ,i1 ,n-start1) (< ,i2 ,n-start2)))
  580.            (declare (fixnum ,i1 ,i2))
  581.            (%bignum-set ,n-dest ,i1
  582.                 (%bignum-ref ,n-src ,i2))))
  583.       `(let ((,n-end1 ,end1)
  584.          (,n-end2 ,end2))
  585.          (do ((,i1 ,start1 (1+ ,i1))
  586.           (,i2 ,start2 (1+ ,i2)))
  587.          ((or (>= ,i1 ,n-end1) (>= ,i2 ,n-end2)))
  588.            (declare (type bignum-index ,i1 ,i2))
  589.            (%bignum-set ,n-dest ,i1
  590.                 (%bignum-ref ,n-src ,i2))))))))
  591.  
  592.  
  593. ;;; WITH-BIGNUM-BUFFERS  --  Internal.
  594. ;;;
  595. ;;; Could do freelisting someday.
  596. ;;;
  597. (defmacro with-bignum-buffers (specs &body body)
  598.   "WITH-BIGNUM-BUFFERS ({(var size [init])}*) Form*"
  599.   (ext:collect ((binds)
  600.         (inits))
  601.     (dolist (spec specs)
  602.       (let ((name (first spec))
  603.         (size (second spec)))
  604.     (binds `(,name (%allocate-bignum ,size)))
  605.     (let ((init (third spec)))
  606.       (when init
  607.         (inits `(bignum-replace ,name ,init))))))
  608.     `(let* ,(binds)
  609.        ,@(inits)
  610.        ,@body)))
  611.  
  612. ) ;EVAL-WHEN
  613.  
  614.  
  615.  
  616. ;;;; GCD.
  617.  
  618. (defun bignum-gcd (a b)
  619.   (declare (type bignum-type a b))
  620.   (let* ((a (if (%bignum-0-or-plusp a (%bignum-length a))
  621.         a
  622.         (negate-bignum a nil)))
  623.      (b (if (%bignum-0-or-plusp b (%bignum-length b))
  624.         b
  625.         (negate-bignum b nil)))
  626.      (len-a (%bignum-length a))
  627.      (len-b (%bignum-length b)))
  628.       (declare (type bignum-index len-a len-b))
  629.     (with-bignum-buffers ((a-buffer len-a a)
  630.               (b-buffer len-b b)
  631.               (res-buffer (max len-a len-b)))
  632.       (let* ((factors-of-two
  633.           (bignum-factors-of-two a-buffer len-a
  634.                      b-buffer len-b))
  635.          (len-a (make-gcd-bignum-odd
  636.              a-buffer
  637.              (bignum-buffer-ashift-right a-buffer len-a
  638.                          factors-of-two)))
  639.          (len-b (make-gcd-bignum-odd
  640.              b-buffer
  641.              (bignum-buffer-ashift-right b-buffer len-b
  642.                          factors-of-two))))
  643.     (declare (type bignum-index len-a len-b))
  644.     (let ((x a-buffer)
  645.           (len-x len-a)
  646.           (y b-buffer)
  647.           (len-y len-b)
  648.           (z res-buffer))
  649.       (loop
  650.         (multiple-value-bind
  651.         (u v len-v r len-r)
  652.         (bignum-gcd-order-and-subtract x len-x y len-y z)
  653.           (declare (type bignum-index len-v len-r))
  654.           (when (and (= len-r 1) (zerop (%bignum-ref r 0)))
  655.         (if (zerop factors-of-two)
  656.             (let ((ret (%allocate-bignum len-v)))
  657.               (dotimes (i len-v)
  658.             (setf (%bignum-ref ret i) (%bignum-ref v i)))
  659.               (return (%normalize-bignum ret len-v)))
  660.             (return (bignum-ashift-left v factors-of-two len-v))))
  661.           (setf x v  len-x len-v)
  662.           (setf y r  len-y (make-gcd-bignum-odd r len-r))
  663.           (setf z u))))))))
  664.  
  665.  
  666. (defun bignum-gcd-order-and-subtract (a len-a b len-b res)
  667.   (declare (type bignum-index len-a len-b) (type bignum-type a b))
  668.   (cond ((= len-a len-b)
  669.      (do ((i (1- len-a) (1- i)))
  670.          ((= i -1)
  671.           (setf (%bignum-ref res 0) 0)
  672.           (values a b len-b res 1))
  673.        (let ((a-digit (%bignum-ref a i))
  674.          (b-digit (%bignum-ref b i)))
  675.          (cond ((%digit-compare a-digit b-digit))
  676.            ((%digit-greater a-digit b-digit)
  677.             (return
  678.              (values a b len-b res
  679.                  (subtract-bignum-buffers a len-a b len-b res))))
  680.            (t
  681.             (return
  682.              (values b a len-a res
  683.                  (subtract-bignum-buffers b len-b a len-a res))))))))
  684.     ((> len-a len-b)
  685.      (values a b len-b res
  686.          (subtract-bignum-buffers a len-a b len-b res)))
  687.     (t
  688.      (values b a len-a res
  689.          (subtract-bignum-buffers b len-b a len-a res)))))
  690.  
  691. (defun make-gcd-bignum-odd (a len-a)
  692.   (declare (type bignum-type a) (type bignum-index len-a))
  693.   (dotimes (index len-a)
  694.     (declare (type bignum-index index))
  695.     (do ((digit (%bignum-ref a index) (%ashr digit 1))
  696.      (increment 0 (1+ increment)))
  697.     ((zerop digit))
  698.       (declare (type (mod 32) increment))
  699.       (when (oddp digit)
  700.     (return-from make-gcd-bignum-odd
  701.              (bignum-buffer-ashift-right a len-a
  702.                          (+ (* index digit-size)
  703.                             increment)))))))
  704.  
  705. (defun bignum-factors-of-two (a len-a b len-b)
  706.   (declare (type bignum-index len-a len-b) (type bignum-type a))
  707.   (do ((i 0 (1+ i))
  708.        (end (min len-a len-b)))
  709.       ((= i end) (error "Unexpected zero bignums?"))
  710.     (declare (type bignum-index i end))
  711.     (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
  712.       (unless (zerop or-digits)
  713.     (return (do ((j 0 (1+ j))
  714.              (or-digits or-digits (%ashr or-digits 1)))
  715.             ((oddp or-digits) (+ (* i digit-size) j))
  716.           (declare (type (mod 32) j))))))))
  717.  
  718.  
  719. ;;;; Negation
  720.  
  721. (eval-when (compile eval)
  722.  
  723. ;;; BIGNUM-NEGATE-LOOP -- Internal.
  724. ;;;
  725. ;;; This negates bignum-len digits of bignum, storing the resulting digits into
  726. ;;; result (possibly EQ to bignum) and returning whatever end-carry there is.
  727. ;;;
  728. (defmacro bignum-negate-loop (bignum bignum-len &optional (result nil resultp))
  729.   (let ((carry (gensym))
  730.     (end (gensym))
  731.     (value (gensym))
  732.     (last (gensym)))
  733.     `(let* (,@(if (not resultp) `(,last))
  734.         (,carry 
  735.          (multiple-value-bind (,value ,carry)
  736.                   (%add-with-carry
  737.                    (%lognot (%bignum-ref ,bignum 0)) 1 0)
  738.            ,(if resultp
  739.             `(setf (%bignum-ref ,result 0) ,value)
  740.             `(setf ,last ,value))
  741.            ,carry))
  742.         (i 1)
  743.         (,end ,bignum-len))
  744.        (declare (type bit ,carry)
  745.         (type bignum-index i ,end))
  746.        (loop
  747.      (when (= i ,end) (return))
  748.      (multiple-value-bind (,value temp)
  749.                   (%add-with-carry
  750.                    (%lognot (%bignum-ref ,bignum i)) 0 ,carry)
  751.        ,(if resultp
  752.         `(setf (%bignum-ref ,result i) ,value)
  753.         `(setf ,last ,value))
  754.        (setf ,carry temp))
  755.      (incf i))
  756.        ,(if resultp carry `(values ,carry ,last)))))
  757.  
  758. ) ;EVAL-WHEN
  759.  
  760. ;;; NEGATE-BIGNUM -- Public.
  761. ;;;
  762. ;;; Fully-normalize is an internal optional.  It cause this to always return
  763. ;;; a bignum, without any extraneous digits, and it never returns a fixnum.
  764. ;;;
  765. (defun negate-bignum (x &optional (fully-normalize t))
  766.   (declare (type bignum-type x))
  767.   (let* ((len-x (%bignum-length x))
  768.      (len-res (1+ len-x))
  769.      (res (%allocate-bignum len-res)))
  770.     (declare (type bignum-index len-x len-res)) ;Test len-res for range?
  771.     (let ((carry (bignum-negate-loop x len-x res)))
  772.       (setf (%bignum-ref res len-x)
  773.         (%add-with-carry (%lognot (%sign-digit x len-x)) 0 carry)))
  774.     (if fully-normalize
  775.     (%normalize-bignum res len-res)
  776.     (%mostly-normalize-bignum res len-res))))
  777.  
  778. ;;; NEGATE-BIGNUM-IN-PLACE -- Internal.
  779. ;;;
  780. ;;; This assumes bignum is positive; that is, the result of negating it will
  781. ;;; stay in the provided allocated bignum.
  782. ;;;
  783. (defun negate-bignum-in-place (bignum)
  784.   (bignum-negate-loop bignum (%bignum-length bignum) bignum)
  785.   bignum)
  786.  
  787.  
  788.  
  789. ;;;; Shifting.
  790.  
  791. (defconstant all-ones-digit #xFFFFFFFF)
  792.  
  793. (eval-when (compile eval)
  794.  
  795. ;;; SHIFT-RIGHT-UNALIGNED -- Internal.
  796. ;;;
  797. ;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
  798. ;;; BIGNUM-LDB-BIGNUM-RES.  They supply a termination form that references
  799. ;;; locals established by this form.  Source is the source bignum.  Start-digit
  800. ;;; is the first digit in source from which we pull bits.  Start-pos is the
  801. ;;; first bit we want.  Res-len-form is the form that computes the length of
  802. ;;; the resulting bignum.  Termination is a DO termination form with a test and
  803. ;;; body.  When result is supplied, it is the variable to which this binds a
  804. ;;; newly allocated bignum.
  805. ;;;
  806. ;;; Given start-pos, 1-31 inclusively, of shift, we form the j'th resulting
  807. ;;; digit from high bits of the i'th source digit and the start-pos number of
  808. ;;; bits from the i+1'th source digit.
  809. ;;;
  810. (defmacro shift-right-unaligned (source start-digit start-pos res-len-form
  811.                  termination
  812.                  &optional result)
  813.   `(let* ((high-bits-in-first-digit (- digit-size ,start-pos))
  814.       (res-len ,res-len-form)
  815.       (res-len-1 (1- res-len))
  816.       ,@(if result `((,result (%allocate-bignum res-len)))))
  817.      (declare (type bignum-index res-len res-len-1))
  818.      (do ((i ,start-digit i+1)
  819.       (i+1 (1+ ,start-digit) (1+ i+1))
  820.       (j 0 (1+ j)))
  821.      ,termination
  822.        (declare (type bignum-index i i+1 j))
  823.        (setf (%bignum-ref ,(if result result source) j)
  824.          (%logior (%digit-logical-shift-right (%bignum-ref ,source i)
  825.                           ,start-pos)
  826.               (%ashl (%bignum-ref ,source i+1)
  827.                  high-bits-in-first-digit))))))
  828.  
  829. ) ;EVAL-WHEN
  830.  
  831.  
  832. ;;; BIGNUM-ASHIFT-RIGHT -- Public.
  833. ;;;
  834. ;;; First compute the number of whole digits to shift, shifting them by
  835. ;;; skipping them when we start to pick up bits, and the number of bits to
  836. ;;; shift the remaining digits into place.  If the number of digits is greater
  837. ;;; than the length of the bignum, then the result is either 0 or -1.  If we
  838. ;;; shift on a digit boundary (that is, n-bits is zero), then we just copy
  839. ;;; digits.  The last branch handles the general case which uses a macro that a
  840. ;;; couple other routines use.  The fifth argument to the macro references
  841. ;;; locals established by the macro.
  842. ;;;
  843. (defun bignum-ashift-right (bignum x)
  844.   (declare (type bignum-type bignum)
  845.        (fixnum x))
  846.   (let ((bignum-len (%bignum-length bignum)))
  847.     (declare (type bignum-index bignum-len))
  848.     (multiple-value-bind (digits n-bits)
  849.              (truncate x digit-size)
  850.       (declare (type bignum-index digits))
  851.       (cond
  852.        ((>= digits bignum-len)
  853.     (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
  854.        ((zerop n-bits)
  855.     (bignum-ashift-right-digits bignum digits))
  856.        (t
  857.     (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
  858.                    ((= j res-len-1)
  859.                 (setf (%bignum-ref res j)
  860.                       (%ashr (%bignum-ref bignum i) n-bits))
  861.                 (%normalize-bignum res res-len))
  862.                    res))))))
  863.  
  864. ;;; BIGNUM-ASHIFT-RIGHT-DIGITS -- Internal.
  865. ;;;
  866. (defun bignum-ashift-right-digits (bignum digits)
  867.   (declare (type bignum-type bignum)
  868.        (type bignum-index digits))
  869.   (let* ((res-len (- (%bignum-length bignum) digits))
  870.      (res (%allocate-bignum res-len)))
  871.     (declare (type bignum-index res-len)
  872.          (type bignum-type res))
  873.     (bignum-replace res bignum :start2 digits)
  874.     (%normalize-bignum res res-len)))
  875.  
  876.  
  877. ;;; BIGNUM-BUFFER-ASHIFT-RIGHT -- Internal.
  878. ;;;
  879. ;;; GCD uses this for an in-place shifting operation.  This is different enough
  880. ;;; from BIGNUM-ASHIFT-RIGHT that it isn't worth folding the bodies into a
  881. ;;; macro, but they share the basic algorithm.  This routine foregoes a first
  882. ;;; test for digits being greater than or equal to bignum-len since that will
  883. ;;; never happen for its uses in GCD.  We did fold the last branch into a macro
  884. ;;; since it was duplicated a few times, and the fifth argument to it
  885. ;;; references locals established by the macro.
  886. ;;;
  887. (defun bignum-buffer-ashift-right (bignum bignum-len x)
  888.   (declare (type bignum-index bignum-len) (fixnum x))
  889.   (multiple-value-bind (digits n-bits)
  890.                (truncate x digit-size)
  891.     (declare (type bignum-index digits))
  892.     (cond
  893.      ((zerop n-bits)
  894.       (let ((new-end (- bignum-len digits)))
  895.     (bignum-replace bignum bignum :end1 new-end :start2 digits
  896.             :end2 bignum-len)
  897.     (%normalize-bignum-buffer bignum new-end)))
  898.      (t
  899.       (shift-right-unaligned bignum digits n-bits (- bignum-len digits)
  900.                  ((= j res-len-1)
  901.                   (setf (%bignum-ref bignum j)
  902.                     (%ashr (%bignum-ref bignum i) n-bits))
  903.                   (%normalize-bignum-buffer bignum res-len)))))))
  904.  
  905.  
  906.  
  907. ;;; BIGNUM-ASHIFT-LEFT -- Public.
  908. ;;;
  909. ;;; This handles shifting a bignum buffer to provide fresh bignum data for some
  910. ;;; internal routines.  We know bignum is safe when called with bignum-len.
  911. ;;; First we compute the number of whole digits to shift, shifting them
  912. ;;; starting to store farther along the result bignum.  If we shift on a digit
  913. ;;; boundary (that is, n-bits is zero), then we just copy digits.  The last
  914. ;;; branch handles the general case.
  915. ;;;
  916. (defun bignum-ashift-left (bignum x &optional bignum-len)
  917.   (declare (type bignum-type bignum)
  918.        (fixnum x)
  919.        (type (or null bignum-index) bignum-len))
  920.   (multiple-value-bind (digits n-bits)
  921.                (truncate x digit-size)
  922.     (let* ((bignum-len (or bignum-len (%bignum-length bignum)))
  923.        (res-len (+ digits bignum-len 1)))
  924.       (when (> res-len maximum-bignum-length)
  925.     (error "Can't represent result of left shift."))
  926.       (if (zerop n-bits)
  927.       (bignum-ashift-left-digits bignum bignum-len digits)
  928.       (bignum-ashift-left-unaligned bignum digits n-bits res-len)))))
  929.  
  930. ;;; BIGNUM-ASHIFT-LEFT-DIGITS -- Internal.
  931. ;;;
  932. (defun bignum-ashift-left-digits (bignum bignum-len digits)
  933.   (declare (type bignum-index bignum-len digits))
  934.   (let* ((res-len (+ bignum-len digits))
  935.      (res (%allocate-bignum res-len)))
  936.     (declare (type bignum-index res-len))
  937.     (bignum-replace res bignum :start1 digits :end1 res-len :end2 bignum-len
  938.             :from-end t)
  939.     res))
  940.  
  941. ;;; BIGNUM-ASHIFT-LEFT-UNALIGNED -- Internal.
  942. ;;;
  943. ;;; BIGNUM-TRUNCATE uses this to store into a bignum buffer by supplying res.
  944. ;;; When res comes in non-nil, then this foregoes allocating a result, and it
  945. ;;; normalizes the buffer instead of the would-be allocated result.
  946. ;;;
  947. ;;; We start storing into one digit higher than digits, storing a whole result
  948. ;;; digit from parts of two contiguous digits from bignum.  When the loop
  949. ;;; finishes, we store the remaining bits from bignum's first digit in the
  950. ;;; first non-zero result digit, digits.  We also grab some left over high
  951. ;;; bits from the last digit of bignum.
  952. ;;; 
  953. (defun bignum-ashift-left-unaligned (bignum digits n-bits res-len
  954.                      &optional (res nil resp))
  955.   (declare (type bignum-index digits res-len)
  956.        (type (mod #.digit-size) n-bits))
  957.   (let* ((remaining-bits (- digit-size n-bits))
  958.      (res-len-1 (1- res-len))
  959.      (res (or res (%allocate-bignum res-len))))
  960.     (declare (type bignum-index res-len res-len-1))
  961.     (do ((i 0 i+1)
  962.      (i+1 1 (1+ i+1))
  963.      (j (1+ digits) (1+ j)))
  964.     ((= j res-len-1)
  965.      (setf (%bignum-ref res digits)
  966.            (%ashl (%bignum-ref bignum 0) n-bits))
  967.      (setf (%bignum-ref res j)
  968.            (%ashr (%bignum-ref bignum i) remaining-bits))
  969.      (if resp
  970.          (%normalize-bignum-buffer res res-len)
  971.          (%normalize-bignum res res-len)))
  972.       (declare (type bignum-index i i+1 j))
  973.       (setf (%bignum-ref res j)
  974.         (%logior (%digit-logical-shift-right (%bignum-ref bignum i)
  975.                          remaining-bits)
  976.              (%ashl (%bignum-ref bignum i+1) n-bits))))))
  977.  
  978.  
  979. ;;;; Relational operators.
  980.  
  981. ;;; BIGNUM-PLUS-P -- Public.
  982. ;;;
  983. ;;; Return T iff bignum is positive.
  984. ;;; 
  985. (defun bignum-plus-p (bignum)
  986.   (declare (type bignum-type bignum))
  987.   (%bignum-0-or-plusp bignum (%bignum-length bignum)))
  988.  
  989. ;;; BIGNUM-COMPARE -- Public.
  990. ;;;
  991. ;;; This compares two bignums returning -1, 0, or 1, depending on whether a
  992. ;;; is less than, equal to, or greater than b.
  993. ;;;
  994. (proclaim '(function bignum-compare (bignum bignum) (integer -1 1)))
  995. (defun bignum-compare (a b)
  996.   (declare (type bignum-type a b))
  997.   (let* ((len-a (%bignum-length a))
  998.      (len-b (%bignum-length b))
  999.      (a-plusp (%bignum-0-or-plusp a len-a))
  1000.      (b-plusp (%bignum-0-or-plusp b len-b)))
  1001.     (declare (type bignum-index len-a len-b))
  1002.     (cond ((not (eq a-plusp b-plusp))
  1003.        (if a-plusp 1 -1))
  1004.       ((= len-a len-b)
  1005.        (do ((i (1- len-a) (1- i)))
  1006.            (())
  1007.          (declare (type bignum-index i))
  1008.          (let ((a-digit (%bignum-ref a i))
  1009.            (b-digit (%bignum-ref b i)))
  1010.            (declare (type bignum-element-type a-digit b-digit))
  1011.            (when (%digit-greater a-digit b-digit)
  1012.          (return 1))
  1013.            (when (%digit-greater b-digit a-digit)
  1014.          (return -1)))
  1015.          (when (zerop i) (return 0))))
  1016.       ((> len-a len-b)
  1017.        (if a-plusp 1 -1))
  1018.       (t (if a-plusp -1 1)))))
  1019.  
  1020.  
  1021. ;;;; Float conversion.
  1022.  
  1023. ;;; FLOAT-BIGNUM-RATIO  --  Internal
  1024. ;;;
  1025. ;;;    Given a Ratio with arbitrarily large numerator and denominator, convert
  1026. ;;; it to a float in the specified Format without causing spurious floating
  1027. ;;; overflows.  What we do is discard all of the numerator and denominator
  1028. ;;; except for the format's precision + guard bits.  We float these
  1029. ;;; integers, do the floating division, and then scaled the result accordingly.
  1030. ;;;
  1031. ;;;    The use of digit-size as the number of guard bits is relatively
  1032. ;;; arbitrary.  We want to keep around some extra bits so that we rarely do
  1033. ;;; round-to-even when there were low bits that could have caused us to round
  1034. ;;; up.  We really only need to discard enough bits to ensure that floating the
  1035. ;;; result doesn't overflow.
  1036. ;;;
  1037. (defun float-bignum-ratio (ratio format)
  1038.   (declare (optimize (ext:inhibit-warnings 3)))
  1039.   (let* ((bits-to-keep (+ (float-format-digits format) digit-size))
  1040.      (num (numerator ratio))
  1041.      (num-len (integer-length num))
  1042.      (num-shift (min (- bits-to-keep num-len) 0))
  1043.      (den (denominator ratio))
  1044.      (den-len (integer-length den))
  1045.      (den-shift (min (- bits-to-keep den-len) 0)))
  1046.     (multiple-value-bind
  1047.     (decoded exp sign)
  1048.     (decode-float (/ (coerce (ash num num-shift) format)
  1049.              (coerce (ash den den-shift) format)))
  1050.      (* sign (scale-float decoded (+ exp (- num-shift) den-shift))))))
  1051.  
  1052.  
  1053. ;;; xxx-FLOAT-FROM-BITS  --  Internal
  1054. ;;;
  1055. ;;;    Make a single or double float with the specified significand, exponent
  1056. ;;; and sign.
  1057. ;;;
  1058. (defun single-float-from-bits (bits exp plusp)
  1059.   (declare (fixnum exp))
  1060.   (declare (optimize (ext:inhibit-warnings 3)))
  1061.   (let ((res (dpb exp
  1062.           vm:single-float-exponent-byte
  1063.           (logandc2 (ext:truly-the (unsigned-byte 31)
  1064.                        (%bignum-ref bits 1))
  1065.                 vm:single-float-hidden-bit))))
  1066.     (make-single-float
  1067.      (if plusp
  1068.      res
  1069.      (logior res (ash -1 vm:float-sign-shift))))))
  1070. ;;;
  1071. (defun double-float-from-bits (bits exp plusp)
  1072.   (declare (fixnum exp))
  1073.   (declare (optimize (ext:inhibit-warnings 3)))
  1074.   (let ((hi (dpb exp
  1075.          vm:double-float-exponent-byte
  1076.          (logandc2 (ext:truly-the (unsigned-byte 31)
  1077.                       (%bignum-ref bits 2))
  1078.                vm:double-float-hidden-bit))))
  1079.     (make-double-float
  1080.      (if plusp
  1081.      hi
  1082.      (logior hi (ash -1 vm:float-sign-shift)))
  1083.      (%bignum-ref bits 1))))
  1084.  
  1085.  
  1086. ;;; BIGNUM-TO-FLOAT   --  Interface
  1087. ;;;
  1088. ;;;    Convert Bignum to a float in the specified Format, rounding to the best
  1089. ;;; approximation.
  1090. ;;;
  1091. (defun bignum-to-float (bignum format)
  1092.   (let* ((plusp (bignum-plus-p bignum))
  1093.      (x (if plusp bignum (negate-bignum bignum)))
  1094.      (len (bignum-integer-length x))
  1095.      (digits (float-format-digits format))
  1096.      (keep (+ digits digit-size))
  1097.      (shift (- keep len))
  1098.      (shifted (if (minusp shift)
  1099.               (bignum-ashift-right x (- shift))
  1100.               (bignum-ashift-left x shift)))
  1101.      (low (%bignum-ref shifted 0))
  1102.      (round-bit (ash 1 (1- digit-size))))
  1103.     (declare (type bignum-index len digits keep) (fixnum shift))
  1104.     (labels ((round-up ()
  1105.            (let ((rounded (add-bignums shifted round-bit)))
  1106.          (if (> (integer-length rounded) keep)
  1107.              (float-from-bits (bignum-ashift-right rounded 1)
  1108.                       (1+ len))
  1109.              (float-from-bits rounded len))))
  1110.          (float-from-bits (bits len)
  1111.            (declare (type bignum-index len))
  1112.            (ecase format
  1113.          (single-float
  1114.           (single-float-from-bits
  1115.            bits
  1116.            (check-exponent len vm:single-float-bias
  1117.                            vm:single-float-normal-exponent-max)
  1118.            plusp))
  1119.          (double-float
  1120.           (double-float-from-bits
  1121.            bits
  1122.            (check-exponent len vm:double-float-bias
  1123.                            vm:double-float-normal-exponent-max)
  1124.            plusp))))
  1125.          (check-exponent (exp bias max)
  1126.            (declare (type bignum-index len))
  1127.            (let ((exp (+ exp bias)))
  1128.          (when (> exp max)
  1129.            (error "Too large to be represented as a ~S:~%  ~S"
  1130.               format x))
  1131.          exp)))
  1132.  
  1133.     (cond
  1134.      ;;
  1135.      ;; Round down if round bit is 0.
  1136.      ((zerop (logand round-bit low))
  1137.       (float-from-bits shifted len))
  1138.      ;;
  1139.      ;; If only round bit is set, then round to even.
  1140.      ((and (= low round-bit)
  1141.        (dotimes (i (- (%bignum-length x) (ceiling keep digit-size))
  1142.                t)
  1143.          (unless (zerop (%bignum-ref x i)) (return nil))))
  1144.       (let ((next (%bignum-ref shifted 1)))
  1145.     (if (oddp next)
  1146.         (round-up)
  1147.         (float-from-bits shifted len))))
  1148.      ;;
  1149.      ;; Otherwise, round up.
  1150.      (t
  1151.       (round-up))))))
  1152.  
  1153.  
  1154. ;;;; Integer length and logcount
  1155.  
  1156. (defun bignum-integer-length (bignum)
  1157.   (declare (type bignum-type bignum))
  1158.   (let* ((len (%bignum-length bignum))
  1159.      (len-1 (1- len))
  1160.      (digit (%bignum-ref bignum len-1)))
  1161.     (declare (type bignum-index len len-1)
  1162.          (type bignum-element-type digit))
  1163.     (+ (integer-length (%fixnum-digit-with-correct-sign digit))
  1164.        (* len-1 digit-size))))
  1165.  
  1166. (defun bignum-logcount (bignum)
  1167.   (declare (type bignum-type bignum))
  1168.   (let* ((length (%bignum-length bignum))
  1169.      (plusp (%bignum-0-or-plusp bignum length))
  1170.      (result 0))
  1171.     (declare (type bignum-index length)
  1172.          (fixnum result))
  1173.     (do ((index 0 (1+ index)))
  1174.     ((= index length) result)
  1175.       (let ((digit (%bignum-ref bignum index)))
  1176.     (declare (type bignum-element-type digit))
  1177.     (incf result (logcount (if plusp digit (%lognot digit))))))))
  1178.  
  1179.  
  1180. ;;;; Logical operations.
  1181.  
  1182. ;;; NOT.
  1183. ;;;
  1184.  
  1185. ;;; BIGNUM-LOGICAL-NOT -- Public.
  1186. ;;;
  1187. (defun bignum-logical-not (a)
  1188.   (declare (type bignum-type a))
  1189.   (let* ((len (%bignum-length a))
  1190.      (res (%allocate-bignum len)))
  1191.     (declare (type bignum-index len))
  1192.     (dotimes (i len res)
  1193.       (declare (type bignum-index i))
  1194.       (setf (%bignum-ref res i) (%lognot (%bignum-ref a i))))))
  1195.  
  1196.  
  1197. ;;; AND.
  1198. ;;;
  1199.  
  1200. ;;; BIGNUM-LOGICAL-AND -- Public.
  1201. ;;;
  1202. (defun bignum-logical-and (a b)
  1203.   (declare (type bignum-type a b))
  1204.   (let* ((len-a (%bignum-length a))
  1205.      (len-b (%bignum-length b))
  1206.      (a-plusp (%bignum-0-or-plusp a len-a))
  1207.      (b-plusp (%bignum-0-or-plusp b len-b)))
  1208.     (declare (type bignum-index len-a len-b))
  1209.     (cond
  1210.      ((< len-a len-b)
  1211.       (if a-plusp
  1212.       (logand-shorter-positive a len-a b (%allocate-bignum len-a))
  1213.       (logand-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
  1214.      ((< len-b len-a)
  1215.       (if b-plusp
  1216.       (logand-shorter-positive b len-b a (%allocate-bignum len-b))
  1217.       (logand-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
  1218.      (t (logand-shorter-positive a len-a b (%allocate-bignum len-a))))))
  1219.  
  1220. ;;; LOGAND-SHORTER-POSITIVE -- Internal.
  1221. ;;;
  1222. ;;; This takes a shorter bignum, a and len-a, that is positive.  Because this
  1223. ;;; is AND, we don't care about any bits longer than a's since its infinite 0
  1224. ;;; sign bits will mask the other bits out of b.  The result is len-a big.
  1225. ;;;
  1226. (defun logand-shorter-positive (a len-a b res)
  1227.   (declare (type bignum-type a b res)
  1228.        (type bignum-index len-a))
  1229.   (dotimes (i len-a)
  1230.     (declare (type bignum-index i))
  1231.     (setf (%bignum-ref res i)
  1232.       (%logand (%bignum-ref a i) (%bignum-ref b i))))
  1233.   (%normalize-bignum res len-a))
  1234.  
  1235. ;;; LOGAND-SHORTER-NEGATIVE -- Internal.
  1236. ;;;
  1237. ;;; This takes a shorter bignum, a and len-a, that is negative.  Because this
  1238. ;;; is AND, we just copy any bits longer than a's since its infinite 1 sign
  1239. ;;; bits will include any bits from b.  The result is len-b big.
  1240. ;;;
  1241. (defun logand-shorter-negative (a len-a b len-b res)
  1242.   (declare (type bignum-type a b res)
  1243.        (type bignum-index len-a len-b))
  1244.   (dotimes (i len-a)
  1245.     (declare (type bignum-index i))
  1246.     (setf (%bignum-ref res i)
  1247.       (%logand (%bignum-ref a i) (%bignum-ref b i))))
  1248.   (do ((i len-a (1+ i)))
  1249.       ((= i len-b))
  1250.     (declare (type bignum-index i))
  1251.     (setf (%bignum-ref res i) (%bignum-ref b i)))
  1252.   (%normalize-bignum res len-b))
  1253.  
  1254. ;;; IOR.
  1255. ;;;
  1256.  
  1257. ;;; BIGNUM-LOGICAL-IOR -- Public.
  1258. ;;;
  1259. (defun bignum-logical-ior (a b)
  1260.   (declare (type bignum-type a b))
  1261.   (let* ((len-a (%bignum-length a))
  1262.      (len-b (%bignum-length b))
  1263.      (a-plusp (%bignum-0-or-plusp a len-a))
  1264.      (b-plusp (%bignum-0-or-plusp b len-b)))
  1265.     (declare (type bignum-index len-a len-b))
  1266.     (cond
  1267.      ((< len-a len-b)
  1268.       (if a-plusp
  1269.       (logior-shorter-positive a len-a b len-b (%allocate-bignum len-b))
  1270.       (logior-shorter-negative a len-a b len-b (%allocate-bignum len-b))))
  1271.      ((< len-b len-a)
  1272.       (if b-plusp
  1273.       (logior-shorter-positive b len-b a len-a (%allocate-bignum len-a))
  1274.       (logior-shorter-negative b len-b a len-a (%allocate-bignum len-a))))
  1275.      (t (logior-shorter-positive a len-a b len-b (%allocate-bignum len-a))))))
  1276.  
  1277. ;;; LOGIOR-SHORTER-POSITIVE -- Internal.
  1278. ;;;
  1279. ;;; This takes a shorter bignum, a and len-a, that is positive.  Because this
  1280. ;;; is IOR, we don't care about any bits longer than a's since its infinite
  1281. ;;; 0 sign bits will mask the other bits out of b out to len-b.  The result
  1282. ;;; is len-b long.
  1283. ;;;
  1284. (defun logior-shorter-positive (a len-a b len-b res)
  1285.   (declare (type bignum-type a b res)
  1286.        (type bignum-index len-a len-b))
  1287.   (dotimes (i len-a)
  1288.     (declare (type bignum-index i))
  1289.     (setf (%bignum-ref res i)
  1290.       (%logior (%bignum-ref a i) (%bignum-ref b i))))
  1291.   (do ((i len-a (1+ i)))
  1292.       ((= i len-b))
  1293.     (declare (type bignum-index i))
  1294.     (setf (%bignum-ref res i) (%bignum-ref b i)))
  1295.   (%normalize-bignum res len-b))
  1296.  
  1297. ;;; LOGIOR-SHORTER-NEGATIVE -- Internal.
  1298. ;;;
  1299. ;;; This takes a shorter bignum, a and len-a, that is negative.  Because this
  1300. ;;; is IOR, we just copy any bits longer than a's since its infinite 1 sign
  1301. ;;; bits will include any bits from b.  The result is len-b long.
  1302. ;;;
  1303. (defun logior-shorter-negative (a len-a b len-b res)
  1304.   (declare (type bignum-type a b res)
  1305.        (type bignum-index len-a len-b))
  1306.   (dotimes (i len-a)
  1307.     (declare (type bignum-index i))
  1308.     (setf (%bignum-ref res i)
  1309.       (%logior (%bignum-ref a i) (%bignum-ref b i))))
  1310.   (do ((i len-a (1+ i))
  1311.        (sign (%sign-digit a len-a)))
  1312.       ((= i len-b))
  1313.     (declare (type bignum-index i))
  1314.     (setf (%bignum-ref res i) sign))
  1315.   (%normalize-bignum res len-b))
  1316.  
  1317. ;;; XOR.
  1318. ;;;
  1319.  
  1320. ;;; BIGNUM-LOGICAL-XOR -- Public.
  1321. ;;;
  1322. (defun bignum-logical-xor (a b)
  1323.   (declare (type bignum-type a b))
  1324.   (let ((len-a (%bignum-length a))
  1325.     (len-b (%bignum-length b)))
  1326.     (declare (type bignum-index len-a len-b))
  1327.     (if (< len-a len-b)
  1328.     (bignum-logical-xor-aux a len-a b len-b (%allocate-bignum len-b))
  1329.     (bignum-logical-xor-aux b len-b a len-a (%allocate-bignum len-a)))))
  1330.  
  1331. ;;; BIGNUM-LOGICAL-XOR-AUX -- Internal.
  1332. ;;;
  1333. ;;; This takes the the shorter of two bignums in a and len-a.  Res is len-b
  1334. ;;; long.  Do the XOR.
  1335. ;;;
  1336. (defun bignum-logical-xor-aux (a len-a b len-b res)
  1337.   (declare (type bignum-type a b res)
  1338.        (type bignum-index len-a len-b))
  1339.   (dotimes (i len-a)
  1340.     (declare (type bignum-index i))
  1341.     (setf (%bignum-ref res i)
  1342.       (%logxor (%bignum-ref a i) (%bignum-ref b i))))
  1343.   (do ((i len-a (1+ i))
  1344.        (sign (%sign-digit a len-a)))
  1345.       ((= i len-b))
  1346.     (declare (type bignum-index i))
  1347.     (setf (%bignum-ref res i) (%logxor sign (%bignum-ref b i))))
  1348.   (%normalize-bignum res len-b))
  1349.  
  1350.  
  1351.  
  1352. ;;;; LDB (load byte)
  1353.  
  1354. #|
  1355. FOR NOW WE DON'T USE LDB OR DPB.  WE USE SHIFTS AND MASKS IN NUMBERS.LISP WHICH
  1356. IS LESS EFFICIENT BUT EASIER TO MAINTAIN.  BILL SAYS THIS CODE CERTAINLY WORKS!
  1357.  
  1358.  
  1359. (defconstant maximum-fixnum-bits #+ibm-rt-pc 27 #-ibm-rt-pc 30)
  1360.  
  1361. ;;; BIGNUM-LOAD-BYTE -- Public.
  1362. ;;;
  1363. (defun bignum-load-byte (byte bignum)
  1364.   (declare (type bignum-type bignum))
  1365.   (let ((byte-len (byte-size byte))
  1366.     (byte-pos (byte-position byte)))
  1367.     (if (< byte-len maximum-fixnum-bits)
  1368.     (bignum-ldb-fixnum-res bignum byte-len byte-pos)
  1369.     (bignum-ldb-bignum-res bignum byte-len byte-pos))))
  1370.  
  1371. ;;; BIGNUM-LDB-FIXNUM-RES -- Internal.
  1372. ;;;
  1373. ;;; This returns a fixnum result of loading a byte from a bignum.  In order, we
  1374. ;;; check for the following conditions:
  1375. ;;;    Insufficient bignum digits to start loading a byte --
  1376. ;;;       Return 0 or byte-len 1's depending on sign of bignum.
  1377. ;;;    One bignum digit containing the whole byte spec --
  1378. ;;;       Grab 'em, shift 'em, and mask out what we don't want.
  1379. ;;;    Insufficient bignum digits to cover crossing a digit boundary --
  1380. ;;;       Grab the available bits in the last digit, and or in whatever
  1381. ;;;       virtual sign bits we need to return a full byte spec.
  1382. ;;;    Else (we cross a digit boundary with all bits available) --
  1383. ;;;       Make a couple masks, grab what we want, shift it around, and
  1384. ;;;       LOGIOR it all together.
  1385. ;;; Because (< maximum-fixnum-bits digit-size) and
  1386. ;;;         (< byte-len maximum-fixnum-bits),
  1387. ;;; we only cross one digit boundary if any.
  1388. ;;;
  1389. (defun bignum-ldb-fixnum-res (bignum byte-len byte-pos)
  1390.   (multiple-value-bind (skipped-digits pos)
  1391.                (truncate byte-pos digit-size)
  1392.     (let ((bignum-len (%bignum-length bignum))
  1393.       (s-digits+1 (1+ skipped-digits)))
  1394.       (declare (type bignum-index bignum-len s-digits+1))
  1395.       (if (>= skipped-digits bignum-len)
  1396.       (if (%bignum-0-or-plusp bignum bignum-len)
  1397.           0
  1398.           (%make-ones byte-len))
  1399.       (let ((end (+ pos byte-len)))
  1400.         (cond ((<= end digit-size)
  1401.            (logand (ash (%bignum-ref bignum skipped-digits) (- pos))
  1402.                ;; Must LOGAND after shift here.
  1403.                (%make-ones byte-len)))
  1404.           ((>= s-digits+1 bignum-len)
  1405.            (let* ((available-bits (- digit-size pos))
  1406.               (res (logand (ash (%bignum-ref bignum skipped-digits)
  1407.                         (- pos))
  1408.                        ;; LOGAND should be unnecessary here
  1409.                        ;; with a logical right shift or a
  1410.                        ;; correct unsigned-byte-32 one.
  1411.                        (%make-ones available-bits))))
  1412.              (if (%bignum-0-or-plusp bignum bignum-len)
  1413.              res
  1414.              (logior (%ashl (%make-ones (- end digit-size))
  1415.                     available-bits)
  1416.                  res))))
  1417.           (t
  1418.            (let* ((high-bits-in-first-digit (- digit-size pos))
  1419.               (high-mask (%make-ones high-bits-in-first-digit))
  1420.               (low-bits-in-next-digit (- end digit-size))
  1421.               (low-mask (%make-ones low-bits-in-next-digit)))
  1422.              (declare (type bignum-element-type high-mask low-mask))
  1423.              (logior (%ashl (logand (%bignum-ref bignum s-digits+1)
  1424.                         low-mask)
  1425.                     high-bits-in-first-digit)
  1426.                  (logand (ash (%bignum-ref bignum skipped-digits)
  1427.                       (- pos))
  1428.                      ;; LOGAND should be unnecessary here with
  1429.                      ;; a logical right shift or a correct
  1430.                      ;; unsigned-byte-32 one.
  1431.                      high-mask))))))))))
  1432.  
  1433. ;;; BIGNUM-LDB-BIGNUM-RES -- Internal.
  1434. ;;;
  1435. ;;; This returns a bignum result of loading a byte from a bignum.  In order, we
  1436. ;;; check for the following conditions:
  1437. ;;;    Insufficient bignum digits to start loading a byte --
  1438. ;;;    Byte-pos starting on a digit boundary --
  1439. ;;;    Byte spec contained in one bignum digit --
  1440. ;;;       Grab the bits we want and stick them in a single digit result.
  1441. ;;;       Since we know byte-pos is non-zero here, we know our single digit
  1442. ;;;       will have a zero high sign bit.
  1443. ;;;    Else (unaligned multiple digits) --
  1444. ;;;       This is like doing a shift right combined with either masking
  1445. ;;;       out unwanted high bits from bignum or filling in virtual sign
  1446. ;;;       bits if bignum had insufficient bits.  We use SHIFT-RIGHT-ALIGNED
  1447. ;;;       and reference lots of local variables this macro establishes.
  1448. ;;;
  1449. (defun bignum-ldb-bignum-res (bignum byte-len byte-pos)
  1450.   (multiple-value-bind (skipped-digits pos)
  1451.                (truncate byte-pos digit-size)
  1452.     (let ((bignum-len (%bignum-length bignum)))
  1453.       (declare (type bignum-index bignum-len))
  1454.       (cond
  1455.        ((>= skipped-digits bignum-len)
  1456.     (make-bignum-virtual-ldb-bits bignum bignum-len byte-len))
  1457.        ((zerop pos)
  1458.     (make-aligned-ldb-bignum bignum bignum-len byte-len skipped-digits))
  1459.        ((< (+ pos byte-len) digit-size)
  1460.     (let ((res (%allocate-bignum 1)))
  1461.       (setf (%bignum-ref res 0)
  1462.         (logand (%ashr (%bignum-ref bignum skipped-digits) pos)
  1463.             (%make-ones byte-len)))
  1464.       res))
  1465.        (t
  1466.     (make-unaligned-ldb-bignum bignum bignum-len
  1467.                    byte-len skipped-digits pos))))))
  1468.  
  1469. ;;; MAKE-BIGNUM-VIRTUAL-LDB-BITS -- Internal.
  1470. ;;;
  1471. ;;; This returns bits from bignum that don't physically exist.  These are
  1472. ;;; all zero or one depending on the sign of the bignum.
  1473. ;;;
  1474. (defun make-bignum-virtual-ldb-bits (bignum bignum-len byte-len)
  1475.   (if (%bignum-0-or-plusp bignum bignum-len)
  1476.       0
  1477.       (multiple-value-bind (res-len-1 extra)
  1478.                (truncate byte-len digit-size)
  1479.     (declare (type bignum-index res-len-1))
  1480.     (let* ((res-len (1+ res-len-1))
  1481.            (res (%allocate-bignum res-len)))
  1482.       (declare (type bignum-index res-len))
  1483.       (do ((j 0 (1+ j)))
  1484.           ((= j res-len-1)
  1485.            (setf (%bignum-ref res j) (%make-ones extra))
  1486.            (%normalize-bignum res res-len))
  1487.         (declare (type bignum-index j))
  1488.         (setf (%bignum-ref res j) all-ones-digit))))))
  1489.  
  1490. ;;; MAKE-ALIGNED-LDB-BIGNUM -- Internal.
  1491. ;;;
  1492. ;;; Since we are picking up aligned digits, we just copy the whole digits
  1493. ;;; we want and fill in extra bits.  We might have a byte-len that extends
  1494. ;;; off the end of the bignum, so we may have to fill in extra 1's if the
  1495. ;;; bignum is negative.
  1496. ;;;
  1497. (defun make-aligned-ldb-bignum (bignum bignum-len byte-len skipped-digits)
  1498.   (multiple-value-bind (res-len-1 extra)
  1499.                (truncate byte-len digit-size)
  1500.     (declare (type bignum-index res-len-1))
  1501.     (let* ((res-len (1+ res-len-1))
  1502.        (res (%allocate-bignum res-len)))
  1503.       (declare (type bignum-index res-len))
  1504.       (do ((i skipped-digits (1+ i))
  1505.        (j 0 (1+ j)))
  1506.       ((or (= j res-len-1) (= i bignum-len))
  1507.        (cond ((< i bignum-len)
  1508.           (setf (%bignum-ref res j)
  1509.             (logand (%bignum-ref bignum i)
  1510.                 (the bignum-element-type (%make-ones extra)))))
  1511.          ((%bignum-0-or-plusp bignum bignum-len))
  1512.          (t
  1513.           (do ((j j (1+ j)))
  1514.               ((= j res-len-1)
  1515.                (setf (%bignum-ref res j) (%make-ones extra)))
  1516.             (setf (%bignum-ref res j) all-ones-digit))))
  1517.        (%normalize-bignum res res-len))
  1518.       (declare (type bignum-index i j))
  1519.       (setf (%bignum-ref res j) (%bignum-ref bignum i))))))
  1520.  
  1521. ;;; MAKE-UNALIGNED-LDB-BIGNUM -- Internal.
  1522. ;;;
  1523. ;;; This grabs unaligned bignum bits from bignum assuming byte-len causes at
  1524. ;;; least one digit boundary crossing.  We use SHIFT-RIGHT-UNALIGNED referencing
  1525. ;;; lots of local variables established by it.
  1526. ;;;
  1527. (defun make-unaligned-ldb-bignum (bignum bignum-len byte-len skipped-digits pos)
  1528.   (multiple-value-bind (res-len-1 extra)
  1529.                (truncate byte-len digit-size)
  1530.     (shift-right-unaligned
  1531.      bignum skipped-digits pos (1+ res-len-1)
  1532.      ((or (= j res-len-1) (= i+1 bignum-len))
  1533.       (cond ((= j res-len-1)
  1534.          (cond
  1535.           ((< extra high-bits-in-first-digit)
  1536.            (setf (%bignum-ref res j)
  1537.              (logand (ash (%bignum-ref bignum i) minus-start-pos)
  1538.                  ;; Must LOGAND after shift here.
  1539.                  (%make-ones extra))))
  1540.           (t
  1541.            (setf (%bignum-ref res j)
  1542.              (logand (ash (%bignum-ref bignum i) minus-start-pos)
  1543.                  ;; LOGAND should be unnecessary here with a logical
  1544.                  ;; right shift or a correct unsigned-byte-32 one.
  1545.                  high-mask))
  1546.            (when (%bignum-0-or-plusp bignum bignum-len)
  1547.          (setf (%bignum-ref res j)
  1548.                (logior (%bignum-ref res j)
  1549.                    (%ashl (%make-ones
  1550.                        (- extra high-bits-in-first-digit))
  1551.                       high-bits-in-first-digit)))))))
  1552.         (t
  1553.          (setf (%bignum-ref res j)
  1554.            (logand (ash (%bignum-ref bignum i) minus-start-pos)
  1555.                ;; LOGAND should be unnecessary here with a logical
  1556.                ;; right shift or a correct unsigned-byte-32 one.
  1557.                high-mask))
  1558.          (unless (%bignum-0-or-plusp bignum bignum-len)
  1559.            ;; Fill in upper half of this result digit with 1's.
  1560.            (setf (%bignum-ref res j)
  1561.              (logior (%bignum-ref res j)
  1562.                  (%ashl low-mask high-bits-in-first-digit)))
  1563.            ;; Fill in any extra 1's we need to be byte-len long.
  1564.            (do ((j (1+ j) (1+ j)))
  1565.            ((>= j res-len-1)
  1566.             (setf (%bignum-ref res j) (%make-ones extra)))
  1567.          (setf (%bignum-ref res j) all-ones-digit)))))
  1568.       (%normalize-bignum res res-len))
  1569.      res)))
  1570.  
  1571.  
  1572.  
  1573. ;;;; DPB (deposit byte).  
  1574.  
  1575. (defun bignum-deposit-byte (new-byte byte-spec bignum)
  1576.   (declare (type bignum-type bignum))
  1577.   (let* ((byte-len (byte-size byte-spec))
  1578.      (byte-pos (byte-position byte-spec))
  1579.      (bignum-len (%bignum-length bignum))
  1580.      (bignum-plusp (%bignum-0-or-plusp bignum bignum-len))
  1581.      (byte-end (+ byte-pos byte-len))
  1582.      (res-len (1+ (max (ceiling byte-end digit-size) bignum-len)))
  1583.      (res (%allocate-bignum res-len)))
  1584.     (declare (type bignum-index bignum-len res-len))
  1585.     ;;
  1586.     ;; Fill in an extra sign digit in case we set what would otherwise be the
  1587.     ;; last digit's last bit.  Normalize at the end in case this was
  1588.     ;; unnecessary.
  1589.     (unless bignum-plusp
  1590.       (setf (%bignum-ref res (1- res-len)) all-ones-digit))
  1591.     (multiple-value-bind (end-digit end-bits)
  1592.              (truncate byte-end digit-size)
  1593.       (declare (type bignum-index end-digit))
  1594.       ;;
  1595.       ;; Fill in bits from bignum up to byte-pos.
  1596.       (multiple-value-bind (pos-digit pos-bits)
  1597.                (truncate byte-pos digit-size)
  1598.     (declare (type bignum-index pos-digit))
  1599.     (do ((i 0 (1+ i))
  1600.          (end (min pos-digit bignum-len)))
  1601.         ((= i end)
  1602.          (cond ((< i bignum-len)
  1603.             (unless (zerop pos-bits)
  1604.               (setf (%bignum-ref res i)
  1605.                 (logand (%bignum-ref bignum i)
  1606.                     (%make-ones pos-bits)))))
  1607.            (bignum-plusp)
  1608.            (t
  1609.             (do ((i i (1+ i)))
  1610.             ((= i pos-digit)
  1611.              (unless (zerop pos-bits)
  1612.                (setf (%bignum-ref res i) (%make-ones pos-bits))))
  1613.               (setf (%bignum-ref res i) all-ones-digit)))))
  1614.       (setf (%bignum-ref res i) (%bignum-ref bignum i)))
  1615.     ;;
  1616.     ;; Fill in bits from new-byte.
  1617.     (if (typep new-byte 'fixnum)
  1618.         (deposit-fixnum-bits new-byte byte-len pos-digit pos-bits
  1619.                  end-digit end-bits res)
  1620.         (deposit-bignum-bits new-byte byte-len pos-digit pos-bits
  1621.                  end-digit end-bits res)))
  1622.       ;;
  1623.       ;; Fill in remaining bits from bignum after byte-spec.
  1624.       (when (< end-digit bignum-len)
  1625.     (setf (%bignum-ref res end-digit)
  1626.           (logior (logand (%bignum-ref bignum end-digit)
  1627.                   (%ashl (%make-ones (- digit-size end-bits))
  1628.                      end-bits))
  1629.               ;; DEPOSIT-FIXNUM-BITS and DEPOSIT-BIGNUM-BITS only store
  1630.               ;; bits from new-byte into res's end-digit element, so
  1631.               ;; we don't need to mask out unwanted high bits.
  1632.               (%bignum-ref res end-digit)))
  1633.     (do ((i (1+ end-digit) (1+ i)))
  1634.         ((= i bignum-len))
  1635.       (setf (%bignum-ref res i) (%bignum-ref bignum i)))))
  1636.     (%normalize-bignum res res-len)))
  1637.  
  1638. ;;; DEPOSIT-FIXNUM-BITS -- Internal.
  1639. ;;;
  1640. ;;; This starts at result's pos-digit skipping pos-bits, and it stores bits
  1641. ;;; from new-byte, a fixnum, into result.  It effectively stores byte-len
  1642. ;;; number of bits, but never stores past end-digit and end-bits in result.
  1643. ;;; The first branch fires when all the bits we want from new-byte are present;
  1644. ;;; if byte-len crosses from the current result digit into the next, the last
  1645. ;;; argument to DEPOSIT-FIXNUM-DIGIT is a mask for those bits.  The second
  1646. ;;; branch handles the need to grab more bits than the fixnum new-byte has, but
  1647. ;;; new-byte is positive; therefore, any virtual bits are zero.  The mask for
  1648. ;;; bits that don't fit in the current result digit is simply the remaining
  1649. ;;; bits in the bignum digit containing new-byte; we don't care if we store
  1650. ;;; some extra in the next result digit since they will be zeros.  The last
  1651. ;;; branch handles the need to grab more bits than the fixnum new-byte has, but
  1652. ;;; new-byte is negative; therefore, any virtual bits must be explicitly filled
  1653. ;;; in as ones.  We call DEPOSIT-FIXNUM-DIGIT to grab what bits actually exist
  1654. ;;; and to fill in the current result digit.
  1655. ;;;
  1656. (defun deposit-fixnum-bits (new-byte byte-len pos-digit pos-bits 
  1657.                 end-digit end-bits result)
  1658.   (declare (type bignum-index pos-digit end-digit))
  1659.   (let ((other-bits (- digit-size pos-bits))
  1660.     (new-byte-digit (%fixnum-to-digit new-byte)))
  1661.     (declare (type bignum-element-type new-byte-digit))
  1662.     (cond ((< byte-len maximum-fixnum-bits)
  1663.        (deposit-fixnum-digit new-byte-digit byte-len pos-digit pos-bits
  1664.                  other-bits result
  1665.                  (- byte-len other-bits)))
  1666.       ((or (plusp new-byte) (zerop new-byte))
  1667.        (deposit-fixnum-digit new-byte-digit byte-len pos-digit pos-bits
  1668.                  other-bits result pos-bits))
  1669.       (t
  1670.        (multiple-value-bind
  1671.            (digit bits)
  1672.            (deposit-fixnum-digit new-byte-digit byte-len pos-digit pos-bits
  1673.                      other-bits result
  1674.                      (if (< (- byte-len other-bits) digit-size)
  1675.                      (- byte-len other-bits)
  1676.                      digit-size))
  1677.          (declare (type bignum-index digit))
  1678.          (cond ((< digit end-digit)
  1679.             (setf (%bignum-ref result digit)
  1680.               (logior (%bignum-ref result digit)
  1681.                   (%ashl (%make-ones (- digit-size bits)) bits)))
  1682.             (do ((i (1+ digit) (1+ i)))
  1683.             ((= i end-digit)
  1684.              (setf (%bignum-ref result i) (%make-ones end-bits)))
  1685.               (setf (%bignum-ref result i) all-ones-digit)))
  1686.            ((> digit end-digit))
  1687.            ((< bits end-bits)
  1688.             (setf (%bignum-ref result digit)
  1689.               (logior (%bignum-ref result digit)
  1690.                   (%ashl (%make-ones (- end-bits bits))
  1691.                      bits))))))))))
  1692.  
  1693. ;;; DEPOSIT-FIXNUM-DIGIT -- Internal.
  1694. ;;;
  1695. ;;; This fills in the current result digit from new-byte-digit.  The first case
  1696. ;;; handles everything we want fitting in the current digit, and other-bits is
  1697. ;;; the number of bits remaining to be filled in result's current digit.  This
  1698. ;;; number is digit-size minus pos-bits.  The second branch handles filling in
  1699. ;;; result's current digit, and it shoves the unused bits of new-byte-digit
  1700. ;;; into the next result digit.  This is correct regardless of new-byte-digit's
  1701. ;;; sign.  It returns the new current result digit and how many bits already
  1702. ;;; filled in the result digit.
  1703. ;;;
  1704. (defun deposit-fixnum-digit (new-byte-digit byte-len pos-digit pos-bits
  1705.                  other-bits result next-digit-bits-needed)
  1706.   (declare (type bignum-index pos-digit)
  1707.        (type bignum-element-type new-byte-digit next-digit-mask))
  1708.   (cond ((<= byte-len other-bits)
  1709.      ;; Bits from new-byte fit in the current result digit.
  1710.      (setf (%bignum-ref result pos-digit)
  1711.            (logior (%bignum-ref result pos-digit)
  1712.                (%ashl (logand new-byte-digit (%make-ones byte-len))
  1713.                   pos-bits)))
  1714.      (if (= byte-len other-bits)
  1715.          (values (1+ pos-digit) 0)
  1716.          (values pos-digit (+ byte-len pos-bits))))
  1717.     (t
  1718.      ;; Some of new-byte's bits go in current result digit.
  1719.      (setf (%bignum-ref result pos-digit)
  1720.            (logior (%bignum-ref result pos-digit)
  1721.                (%ashl (logand new-byte-digit (%make-ones other-bits))
  1722.                   pos-bits)))
  1723.      (let ((pos-digit+1 (1+ pos-digit)))
  1724.        ;; The rest of new-byte's bits go in the next result digit.
  1725.        (setf (%bignum-ref result pos-digit+1)
  1726.          (logand (ash new-byte-digit (- other-bits))
  1727.              ;; Must LOGAND after shift here.
  1728.              (%make-ones next-digit-bits-needed)))
  1729.        (if (= next-digit-bits-needed digit-size)
  1730.            (values (1+ pos-digit+1) 0)
  1731.            (values pos-digit+1 next-digit-bits-needed))))))
  1732.  
  1733. ;;; DEPOSIT-BIGNUM-BITS -- Internal.
  1734. ;;;
  1735. ;;; This starts at result's pos-digit skipping pos-bits, and it stores bits
  1736. ;;; from new-byte, a bignum, into result.  It effectively stores byte-len
  1737. ;;; number of bits, but never stores past end-digit and end-bits in result.
  1738. ;;; When handling a starting bit unaligned with a digit boundary, we check
  1739. ;;; in the second branch for the byte spec fitting into the pos-digit element
  1740. ;;; after after pos-bits; DEPOSIT-UNALIGNED-BIGNUM-BITS expects at least one
  1741. ;;; digit boundary crossing.
  1742. ;;;
  1743. (defun deposit-bignum-bits (bignum-byte byte-len pos-digit pos-bits 
  1744.                 end-digit end-bits result)
  1745.   (declare (type bignum-index pos-digit end-digit))
  1746.   (cond ((zerop pos-bits)
  1747.      (deposit-aligned-bignum-bits bignum-byte pos-digit end-digit end-bits
  1748.                       result))
  1749.     ((or (= end-digit pos-digit)
  1750.          (and (= end-digit (1+ pos-digit))
  1751.           (zerop end-bits)))
  1752.      (setf (%bignum-ref result pos-digit)
  1753.            (logior (%bignum-ref result pos-digit)
  1754.                (%ashl (logand (%bignum-ref bignum-byte 0)
  1755.                       (%make-ones byte-len))
  1756.                   pos-bits))))
  1757.     (t (deposit-unaligned-bignum-bits bignum-byte pos-digit pos-bits
  1758.                       end-digit end-bits result))))
  1759.  
  1760. ;;; DEPOSIT-ALIGNED-BIGNUM-BITS -- Internal.
  1761. ;;;
  1762. ;;; This deposits bits from bignum-byte into result starting at pos-digit and
  1763. ;;; the zero'th bit.  It effectively only stores bits to end-bits in the
  1764. ;;; end-digit element of result.  The loop termination code takes care of
  1765. ;;; picking up the last digit's bits or filling in virtual negative sign bits.
  1766. ;;;
  1767. (defun deposit-aligned-bignum-bits (bignum-byte pos-digit end-digit end-bits
  1768.                     result)
  1769.   (declare (type bignum-index pos-digit end-digit))
  1770.   (let* ((bignum-len (%bignum-length bignum-byte))
  1771.      (bignum-plusp (%bignum-0-or-plusp bignum-byte bignum-len)))
  1772.     (declare (type bignum-index bignum-len))
  1773.     (do ((i 0 (1+ i ))
  1774.      (j pos-digit (1+ j)))
  1775.     ((or (= j end-digit) (= i bignum-len))
  1776.      (cond ((= j end-digit)
  1777.         (cond ((< i bignum-len)
  1778.                (setf (%bignum-ref result j)
  1779.                  (logand (%bignum-ref bignum-byte i)
  1780.                      (%make-ones end-bits))))
  1781.               (bignum-plusp)
  1782.               (t
  1783.                (setf (%bignum-ref result j) (%make-ones end-bits)))))
  1784.            (bignum-plusp)
  1785.            (t
  1786.         (do ((j j (1+ j)))
  1787.             ((= j end-digit)
  1788.              (setf (%bignum-ref result j) (%make-ones end-bits)))
  1789.           (setf (%bignum-ref result j) all-ones-digit)))))
  1790.       (setf (%bignum-ref result j) (%bignum-ref bignum-byte i)))))
  1791.  
  1792. ;;; DEPOSIT-UNALIGNED-BIGNUM-BITS -- Internal.
  1793. ;;;
  1794. ;;; This assumes at least one digit crossing.
  1795. ;;;
  1796. (defun deposit-unaligned-bignum-bits (bignum-byte pos-digit pos-bits
  1797.                       end-digit end-bits result)
  1798.   (declare (type bignum-index pos-digit end-digit))
  1799.   (let* ((bignum-len (%bignum-length bignum-byte))
  1800.      (bignum-plusp (%bignum-0-or-plusp bignum-byte bignum-len))
  1801.      (low-mask (%make-ones pos-bits))
  1802.      (bits-past-pos-bits (- digit-size pos-bits))
  1803.      (high-mask (%make-ones bits-past-pos-bits))
  1804.      (minus-high-bits (- bits-past-pos-bits)))
  1805.     (declare (type bignum-element-type low-mask high-mask)
  1806.          (type bignum-index bignum-len))
  1807.     (do ((i 0 (1+ i))
  1808.      (j pos-digit j+1)
  1809.      (j+1 (1+ pos-digit) (1+ j+1)))
  1810.     ((or (= j end-digit) (= i bignum-len))
  1811.      (cond
  1812.       ((= j end-digit)
  1813.        (setf (%bignum-ref result j)
  1814.          (cond
  1815.           ((>= pos-bits end-bits)
  1816.            (logand (%bignum-ref result j) (%make-ones end-bits)))
  1817.           ((< i bignum-len)
  1818.            (logior (%bignum-ref result j)
  1819.                (%ashl (logand (%bignum-ref bignum-byte i)
  1820.                       (%make-ones (- end-bits pos-bits)))
  1821.                   pos-bits)))
  1822.           (bignum-plusp
  1823.            (logand (%bignum-ref result j)
  1824.                ;; 0's between pos-bits and end-bits positions.
  1825.                (logior (%ashl (%make-ones (- digit-size end-bits))
  1826.                       end-bits)
  1827.                    low-mask)))
  1828.           (t (logior (%bignum-ref result j)
  1829.                  (%ashl (%make-ones (- end-bits pos-bits))
  1830.                     pos-bits))))))
  1831.       (bignum-plusp)
  1832.       (t
  1833.        (setf (%bignum-ref result j)
  1834.          (%ashl (%make-ones bits-past-pos-bits) pos-bits))
  1835.        (do ((j j+1 (1+ j)))
  1836.            ((= j end-digit)
  1837.         (setf (%bignum-ref result j) (%make-ones end-bits)))
  1838.          (declare (type bignum-index j))
  1839.          (setf (%bignum-ref result j) all-ones-digit)))))
  1840.       (declare (type bignum-index i j j+1))
  1841.       (let ((digit (%bignum-ref bignum-byte i)))
  1842.     (declare (type bignum-element-type digit))
  1843.     (setf (%bignum-ref result j)
  1844.           (logior (%bignum-ref result j)
  1845.               (%ashl (logand digit high-mask) pos-bits)))
  1846.     (setf (%bignum-ref result j+1)
  1847.           (logand (ash digit minus-high-bits)
  1848.               ;; LOGAND should be unnecessary here with a logical right
  1849.               ;; shift or a correct unsigned-byte-32 one.
  1850.               low-mask))))))
  1851.  
  1852. |#
  1853.  
  1854.  
  1855. ;;;; TRUNCATE.
  1856.  
  1857. ;;; This is the original sketch of the algorithm from which I implemented this
  1858. ;;; TRUNCATE, assuming both operands are bignums.  I should modify this to work
  1859. ;;; with the documentation on my functions, as a general introduction.  I've
  1860. ;;; left this here just in case someone needs it in the future.  Don't look
  1861. ;;; at this unless reading the functions' comments leaves you at a loss.
  1862. ;;; Remember this comes from Knuth, so the book might give you the right general
  1863. ;;; overview.
  1864. ;;; 
  1865. ;;;
  1866. ;;; (truncate x y):
  1867. ;;;
  1868. ;;; If X's magnitude is less than Y's, then result is 0 with remainder X.
  1869. ;;;
  1870. ;;; Make x and y positive, copying x if it is already positive.
  1871. ;;;
  1872. ;;; Shift y left until there's a 1 in the 30'th bit (most significant, non-sign
  1873. ;;;       digit)
  1874. ;;;    Just do most sig digit to determine how much to shift whole number.
  1875. ;;; Shift x this much too.
  1876. ;;; Remember this initial shift count.
  1877. ;;;
  1878. ;;; Allocate q to be len-x minus len-y quantity plus 1.
  1879. ;;;
  1880. ;;; i = last digit of x.
  1881. ;;; k = last digit of q.
  1882. ;;;
  1883. ;;; LOOP
  1884. ;;;
  1885. ;;; j = last digit of y.
  1886. ;;;
  1887. ;;; compute guess.
  1888. ;;; if x[i] = y[j] then g = #xFFFFFFFF
  1889. ;;; else g = x[i]x[i-1]/y[j].
  1890. ;;;
  1891. ;;; check guess.
  1892. ;;; %UNSIGNED-MULTIPLY returns b and c defined below.
  1893. ;;;    a = x[i-1] - (logand (* g y[j]) #xFFFFFFFF).
  1894. ;;;       Use %UNSIGNED-MULTIPLY taking low-order result.
  1895. ;;;    b = (logand (ash (* g y[j-1]) -32) #xFFFFFFFF).
  1896. ;;;    c = (logand (* g y[j-1]) #xFFFFFFFF).
  1897. ;;; if a < b, okay.
  1898. ;;; if a > b, guess is too high
  1899. ;;;    g = g - 1; go back to "check guess".
  1900. ;;; if a = b and c > x[i-2], guess is too high
  1901. ;;;    g = g - 1; go back to "check guess".
  1902. ;;; GUESS IS 32-BIT NUMBER, SO USE THING TO KEEP IN SPECIAL REGISTER
  1903. ;;; SAME FOR A, B, AND C.
  1904. ;;;
  1905. ;;; Subtract g * y from x[i - len-y+1]..x[i].  See paper for doing this in step.
  1906. ;;; If x[i] < 0, guess is fucked.
  1907. ;;;    negative g, then add 1
  1908. ;;;    zero or positive g, then subtract 1
  1909. ;;; AND add y back into x[len-y+1..i].
  1910. ;;;
  1911. ;;; q[k] = g.
  1912. ;;; i = i - 1.
  1913. ;;; k = k - 1.
  1914. ;;;
  1915. ;;; If k>=0, goto LOOP.
  1916. ;;;
  1917. ;;;
  1918. ;;; Now quotient is good, but remainder is not.
  1919. ;;; Shift x right by saved initial left shifting count.
  1920. ;;;
  1921. ;;; Check quotient and remainder signs.
  1922. ;;; x pos y pos --> q pos r pos
  1923. ;;; x pos y neg --> q neg r pos
  1924. ;;; x neg y pos --> q neg r neg
  1925. ;;; x neg y neg --> q pos r neg
  1926. ;;;
  1927. ;;; Normalize quotient and remainder.  Cons result if necessary.
  1928. ;;;
  1929.  
  1930.  
  1931.  
  1932. ;;; These are used by BIGNUM-TRUNCATE and friends in the general case.
  1933. ;;;
  1934. (defvar *truncate-x*)
  1935. (defvar *truncate-y*)
  1936.  
  1937. ;;; BIGNUM-TRUNCATE -- Public.
  1938. ;;;
  1939. ;;; This divides x by y returning the quotient and remainder.  In the general
  1940. ;;; case, we shift y to setup for the algorithm, and we use two buffers to save
  1941. ;;; consing intermediate values.  X gets destructively modified to become the
  1942. ;;; remainder, and we have to shift it to account for the initial Y shift.
  1943. ;;; After we multiple bind q and r, we first fix up the signs and then return
  1944. ;;; the normalized results.
  1945. ;;;
  1946. (defun bignum-truncate (x y)
  1947.   (declare (type bignum-type x y))
  1948.   (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
  1949.      (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
  1950.      (x (if x-plusp x (negate-bignum x nil)))
  1951.      (y (if y-plusp y (negate-bignum y nil)))
  1952.      (len-x (%bignum-length x))
  1953.      (len-y (%bignum-length y)))
  1954.     (multiple-value-bind
  1955.     (q r)
  1956.     (cond ((< len-y 2)
  1957.            (bignum-truncate-single-digit x len-x y))
  1958.           ((plusp (bignum-compare y x))
  1959.            (let ((res (%allocate-bignum len-x)))
  1960.          (dotimes (i len-x)
  1961.            (setf (%bignum-ref res i) (%bignum-ref x i)))
  1962.          (values 0 res)))
  1963.           (t
  1964.            (let ((len-x+1 (1+ len-x)))
  1965.          (with-bignum-buffers ((*truncate-x* len-x+1)
  1966.                        (*truncate-y* (1+ len-y)))
  1967.            (let ((y-shift (shift-y-for-truncate y)))
  1968.              (shift-and-store-truncate-buffers x len-x y len-y y-shift)
  1969.              (values (do-truncate len-x+1 len-y)
  1970.                  ;; DO-TRUNCATE must execute first.
  1971.                  (cond
  1972.                   ((zerop y-shift)
  1973.                    (let ((res (%allocate-bignum len-y)))
  1974.                  (declare (type bignum-type res))
  1975.                  (bignum-replace res *truncate-x* :end2 len-y)
  1976.                  (%normalize-bignum res len-y)))
  1977.                   (t
  1978.                    (shift-right-unaligned
  1979.                 *truncate-x* 0 y-shift len-y
  1980.                 ((= j res-len-1)
  1981.                  (setf (%bignum-ref res j)
  1982.                        (%ashr (%bignum-ref *truncate-x* i)
  1983.                           y-shift))
  1984.                  (%normalize-bignum res res-len))
  1985.                 res)))))))))
  1986.       (let ((quotient (cond ((eq x-plusp y-plusp) q)
  1987.                 ((typep q 'fixnum) (the fixnum (- q)))
  1988.                 (t (negate-bignum-in-place q))))
  1989.         (rem (cond (x-plusp r)
  1990.                ((typep r 'fixnum) (the fixnum (- r)))
  1991.                (t (negate-bignum-in-place r)))))
  1992.     (values (if (typep quotient 'fixnum)
  1993.             quotient
  1994.             (%normalize-bignum quotient (%bignum-length quotient)))
  1995.         (if (typep rem 'fixnum)
  1996.             rem
  1997.             (%normalize-bignum rem (%bignum-length rem))))))))
  1998.  
  1999. ;;; BIGNUM-TRUNCATE-SINGLE-DIGIT -- Internal.
  2000. ;;;
  2001. ;;; This divides x by y when y is a single bignum digit.  BIGNUM-TRUNCATE fixes
  2002. ;;; up the quotient and remainder with respect to sign and normalization.
  2003. ;;;
  2004. ;;; We don't have to worry about shifting y to make its most significant digit
  2005. ;;; sufficiently large for %FLOOR to return 32-bit quantities for the q-digit
  2006. ;;; and r-digit.  If y is a single digit bignum, it is already large enough
  2007. ;;; for %FLOOR.  That is, it has some bits on pretty high in the digit.
  2008. ;;;
  2009. (defun bignum-truncate-single-digit (x len-x y)
  2010.   (declare (type bignum-index len-x))
  2011.   (let ((q (%allocate-bignum len-x))
  2012.     (r 0)
  2013.     (y (%bignum-ref y 0)))
  2014.     (declare (type bignum-element-type r y))
  2015.     (do ((i (1- len-x) (1- i)))
  2016.     ((minusp i))
  2017.       (multiple-value-bind (q-digit r-digit)
  2018.                (%floor r (%bignum-ref x i) y)
  2019.     (declare (type bignum-element-type q-digit r-digit))
  2020.     (setf (%bignum-ref q i) q-digit)
  2021.     (setf r r-digit)))
  2022.     (let ((rem (%allocate-bignum 1)))
  2023.       (setf (%bignum-ref rem 0) r)
  2024.       (values q rem))))
  2025.  
  2026. ;;; DO-TRUNCATE -- Internal.
  2027. ;;;
  2028. ;;; This divides *truncate-x* by *truncate-y*, and len-x and len-y tell us how
  2029. ;;; much of the buffers we care about.  TRY-BIGNUM-TRUNCATE-GUESS modifies
  2030. ;;; *truncate-x* on each interation, and this buffer becomes our remainder.
  2031. ;;;
  2032. ;;; *truncate-x* definitely has at least three digits, and it has one more than
  2033. ;;; *truncate-y*.  This keeps i, i-1, i-2, and low-x-digit happy.  Thanks to
  2034. ;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS.
  2035. ;;;
  2036. (defun do-truncate (len-x len-y)
  2037.   (declare (type bignum-index len-x len-y))
  2038.   (let* ((len-q (- len-x len-y))
  2039.      ;; Add one for extra sign digit in case high bit is on.
  2040.      (q (%allocate-bignum (1+ len-q)))
  2041.      (k (1- len-q))
  2042.      (y1 (%bignum-ref *truncate-y* (1- len-y)))
  2043.      (y2 (%bignum-ref *truncate-y* (- len-y 2)))
  2044.      (i (1- len-x))
  2045.      (i-1 (1- i))
  2046.      (i-2 (1- i-1))
  2047.      (low-x-digit (- i len-y)))
  2048.     (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
  2049.          (type bignum-element-type y1 y2))
  2050.     (loop
  2051.       (setf (%bignum-ref q k)
  2052.         (try-bignum-truncate-guess
  2053.          ;; This modifies *truncate-x*.  Must access elements each pass.
  2054.          (bignum-truncate-guess y1 y2
  2055.                     (%bignum-ref *truncate-x* i)
  2056.                     (%bignum-ref *truncate-x* i-1)
  2057.                     (%bignum-ref *truncate-x* i-2))
  2058.          len-y low-x-digit))
  2059.       (cond ((zerop k) (return))
  2060.         (t (decf k)
  2061.            (decf low-x-digit)
  2062.            (shiftf i i-1 i-2 (1- i-2)))))
  2063.     q))
  2064.  
  2065. ;;; TRY-BIGNUM-TRUNCATE-GUESS -- Internal.
  2066. ;;;
  2067. ;;; This takes a digit guess, multiplies it by *truncate-y* for a result one
  2068. ;;; greater in length than len-y, and subtracts this result from *truncate-x*.
  2069. ;;; Low-x-digit is the first digit of x to start the subtraction, and we know x
  2070. ;;; is long enough to subtract a len-y plus one length bignum from it.  Next we
  2071. ;;; check the result of the subtraction, and if the high digit in x became
  2072. ;;; negative, then our guess was one too big.  In this case, return one less
  2073. ;;; than guess passed in, and add one value of y back into x to account for
  2074. ;;; subtracting one too many.  Knuth shows that the guess is wrong on the order
  2075. ;;; of 3/b, where b is the base (2 to the digit-size power) -- pretty rarely.
  2076. ;;;
  2077. (defun try-bignum-truncate-guess (guess len-y low-x-digit)
  2078.   (declare (type bignum-index low-x-digit len-y)
  2079.        (type bignum-element-type guess))
  2080.   (let ((carry-digit 0)
  2081.     (borrow 1)
  2082.     (i low-x-digit))
  2083.     (declare (type bignum-element-type carry-digit)
  2084.          (type bignum-index i)
  2085.          (fixnum borrow))
  2086.     ;; Multiply guess and divisor, subtracting from dividend simultaneously.
  2087.     (dotimes (j len-y)
  2088.       (multiple-value-bind (high-digit low-digit)
  2089.                (%multiply-and-add guess (%bignum-ref *truncate-y* j)
  2090.                           carry-digit)
  2091.     (declare (type bignum-element-type high-digit low-digit))
  2092.     (setf carry-digit high-digit)
  2093.     (multiple-value-bind (x temp-borrow)
  2094.                  (%subtract-with-borrow (%bignum-ref *truncate-x* i)
  2095.                             low-digit borrow)
  2096.       (declare (type bignum-element-type x)
  2097.            (fixnum temp-borrow))
  2098.       (setf (%bignum-ref *truncate-x* i) x)
  2099.       (setf borrow temp-borrow)))
  2100.       (incf i))
  2101.     (setf (%bignum-ref *truncate-x* i)
  2102.       (%subtract-with-borrow (%bignum-ref *truncate-x* i)
  2103.                  carry-digit borrow))
  2104.     ;; See if guess is off by one, adding one Y back in if necessary.
  2105.     (cond ((%digit-0-or-plusp (%bignum-ref *truncate-x* i))
  2106.        guess)
  2107.       (t
  2108.        ;; If subtraction has negative result, add one divisor value back
  2109.        ;; in.  The guess was one too large in magnitude.
  2110.        (let ((i low-x-digit)
  2111.          (carry 0))
  2112.          (dotimes (j len-y)
  2113.            (multiple-value-bind (v k)
  2114.                     (%add-with-carry (%bignum-ref *truncate-y* j)
  2115.                              (%bignum-ref *truncate-x* i)
  2116.                              carry)
  2117.          (declare (type bignum-element-type v))
  2118.          (setf (%bignum-ref *truncate-x* i) v)
  2119.          (setf carry k))
  2120.            (incf i))
  2121.          (setf (%bignum-ref *truncate-x* i)
  2122.            (%add-with-carry (%bignum-ref *truncate-x* i) 0 carry)))
  2123.        (%subtract-with-borrow guess 1 1)))))
  2124.  
  2125. ;;; BIGNUM-TRUNCATE-GUESS -- Internal.
  2126. ;;;
  2127. ;;; This returns a guess for the next division step.  Y1 is the highest y
  2128. ;;; digit, and y2 is the second to highest y digit.  The x... variables are
  2129. ;;; the three highest x digits for the next division step.
  2130. ;;;
  2131. ;;; From Knuth, our guess is either all ones or x-i and x-i-1 divided by y1,
  2132. ;;; depending on whether x-i and y1 are the same.  We test this guess by
  2133. ;;; determining whether guess*y2 is greater than the three high digits of x
  2134. ;;; minus guess*y1 shifted left one digit:
  2135. ;;;    ------------------------------
  2136. ;;;   |    x-i    |   x-i-1  | x-i-2 |
  2137. ;;;    ------------------------------
  2138. ;;;    ------------------------------
  2139. ;;; - | g*y1 high | g*y1 low |   0   |
  2140. ;;;    ------------------------------
  2141. ;;;                ...                   <   guess*y2     ???
  2142. ;;; If guess*y2 is greater, then we decrement our guess by one and try again.
  2143. ;;; This returns a guess that is either correct or one too large.
  2144. ;;;
  2145. (defun bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
  2146.   (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
  2147.   (let ((guess (if (%digit-compare x-i y1)
  2148.            all-ones-digit
  2149.            (%floor x-i x-i-1 y1))))
  2150.     (declare (type bignum-element-type guess))
  2151.     (loop
  2152.       (multiple-value-bind (high-guess*y1 low-guess*y1)
  2153.                (%multiply guess y1)
  2154.     (declare (type bignum-element-type low-guess*y1 high-guess*y1))
  2155.     (multiple-value-bind (high-guess*y2 low-guess*y2)
  2156.                  (%multiply guess y2)
  2157.       (declare (type bignum-element-type high-guess*y2 low-guess*y2))
  2158.       (multiple-value-bind (middle-digit borrow)
  2159.                    (%subtract-with-borrow x-i-1 low-guess*y1 1)
  2160.         (declare (type bignum-element-type middle-digit)
  2161.              (fixnum borrow))
  2162.         ;; Supplying borrow of 1 means there was no borrow, and we know
  2163.         ;; x-i-2 minus 0 requires no borrow.
  2164.         (let ((high-digit (%subtract-with-borrow x-i high-guess*y1 borrow)))
  2165.           (declare (type bignum-element-type high-digit))
  2166.           (if (and (%digit-compare high-digit 0)
  2167.                (or (%digit-greater high-guess*y2 middle-digit)
  2168.                (and (%digit-compare middle-digit high-guess*y2)
  2169.                 (%digit-greater low-guess*y2 x-i-2))))
  2170.           (setf guess (%subtract-with-borrow guess 1 1))
  2171.           (return guess)))))))))
  2172.  
  2173. ;;; SHIFT-Y-FOR-TRUNCATE -- Internal.
  2174. ;;;
  2175. ;;; This returns the amount to shift y to place a one in the second highest
  2176. ;;; bit.  Y must be positive.  If the last digit of y is zero, then y has a
  2177. ;;; one in the previous digit's sign bit, so we know it will take one less
  2178. ;;; than digit-size to get a one where we want.  Otherwise, we count how many
  2179. ;;; right shifts it takes to get zero; subtracting this value from digit-size
  2180. ;;; tells us how many high zeros there are which is one more than the shift
  2181. ;;; amount sought.
  2182. ;;;
  2183. ;;; Note: This is exactly the same as one less than the integer-length of the
  2184. ;;; last digit subtracted from the digit-size.
  2185. ;;; 
  2186. ;;; We shift y to make it sufficiently large that doing the 64-bit by 32-bit
  2187. ;;; %FLOOR calls ensures the quotient and remainder fit in 32-bits.
  2188. ;;;
  2189. (defun shift-y-for-truncate (y)
  2190.   (let* ((len (%bignum-length y))
  2191.      (last (%bignum-ref y (1- len))))
  2192.     (declare (type bignum-index len)
  2193.          (type bignum-element-type last))
  2194.     (- digit-size (integer-length last) 1)))
  2195.  
  2196. ;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS -- Internal.
  2197. ;;;
  2198. ;;; Stores two bignums into the truncation bignum buffers, shifting them on the
  2199. ;;; way in.  This assumes x and y are positive and at least two in length, and
  2200. ;;; it assumes *truncate-x* and *truncate-y* are one digit longer than x and y.
  2201. ;;;
  2202. (defun shift-and-store-truncate-buffers (x len-x y len-y shift)
  2203.   (declare (type bignum-index len-x len-y)
  2204.        (type (integer 0 (#.digit-size)) shift))
  2205.   (cond ((zerop shift)
  2206.      (bignum-replace *truncate-x* x :end1 len-x)
  2207.      (bignum-replace *truncate-y* y :end1 len-y))
  2208.     (t
  2209.      (bignum-ashift-left-unaligned x 0 shift (1+ len-x) *truncate-x*)
  2210.      (bignum-ashift-left-unaligned y 0 shift (1+ len-y) *truncate-y*))))
  2211.  
  2212.  
  2213.  
  2214. ;;;; %FLOOR primitive for BIGNUM-TRUNCATE.
  2215.  
  2216. ;;; When a machine leaves out a 64-bit by 32-bit divide instruction (that is,
  2217. ;;; two bignum-digits divided by one), we have to roll our own (the hard way).
  2218. ;;; Basically, we treat the operation as four 16-bit digits divided by two
  2219. ;;; 16-bit digits.  This means we have duplicated most of the code above to do
  2220. ;;; this nearly general 16-bit digit bignum divide, but we've unrolled loops
  2221. ;;; and made use of other properties of this specific divide situation.
  2222. ;;;
  2223.  
  2224.  
  2225. ;;;
  2226. ;;; %FLOOR for machines with a 32x32 divider.
  2227. ;;;
  2228.  
  2229. (proclaim '(inline 32x16-subtract-with-borrow 32x16-add-with-carry
  2230.            32x16-divide 32x16-multiply 32x16-multiply-split))
  2231.  
  2232. #+32x16-divide
  2233. (defconstant 32x16-base-1 #xFFFF)
  2234.  
  2235. ;;; 32X16-SUBTRACT-WITH-BORROW -- Internal.    [optionally IN ASSEMBLER]
  2236. ;;;
  2237. ;;; This is similar to %SUBTRACT-WITH-BORROW.  It returns a 16-bit difference
  2238. ;;; and a borrow.  Returning a 1 for the borrow means there was no borrow, and
  2239. ;;; 0 means there was one.
  2240. ;;;
  2241. #+32x16-divide
  2242. (defun 32x16-subtract-with-borrow (a b borrow)
  2243.   (declare (type (unsigned-byte 16) a b)
  2244.        (type (integer 0 1) borrow))
  2245.   (let ((diff (+ (- a b) borrow 32x16-base-1)))
  2246.     (declare (type (unsigned-byte 17) diff))
  2247.     (values (logand diff #xFFFF)
  2248.         (ash diff -16))))
  2249.  
  2250. ;;; 32X16-ADD-WITH-CARRY -- Internal.        [optionally IN ASSEMBLER]
  2251. ;;;
  2252. ;;; This adds a and b, 16-bit quantities, with the carry k.  It returns a
  2253. ;;; 16-bit sum and a second value, 0 or 1, indicating whether there was a
  2254. ;;; carry.
  2255. ;;;
  2256. #+32x16-divide
  2257. (defun 32x16-add-with-carry (a b k)
  2258.   (declare (type (unsigned-byte 16) a b)
  2259.        (type (integer 0 1) k))
  2260.   (let ((res (the fixnum (+ a b k))))
  2261.     (declare (type (unsigned-byte 17) res))
  2262.     (if (zerop (the fixnum (logand #x10000 res)))
  2263.     (values res 0)
  2264.     (values (the (unsigned-byte 16) (logand #xFFFF res))
  2265.         1))))
  2266.  
  2267. ;;; 32x16-DIVIDE  --  Internal        [IN ASSEMBLER]
  2268. ;;;
  2269. ;;; This is probably a 32-bit by 32-bit divide instruction.
  2270. ;;;
  2271. #+32x16-divide
  2272. (defun 32x16-divide (a b c)
  2273.   (declare (type (unsigned-byte 16) a b c))
  2274.   (floor (the bignum-element-type
  2275.           (logior (the bignum-element-type (ash a 16))
  2276.               b))
  2277.      c))
  2278.  
  2279. ;;; 32X16-MULTIPLY -- Internal.        [optionally IN ASSEMBLER]
  2280. ;;;
  2281. ;;; This basically exists since we know the answer won't overflow
  2282. ;;; bignum-element-type.  It's probably just a basic multiply instruction, but
  2283. ;;; it can't cons an intermediate bignum.  The result goes in a non-descriptor
  2284. ;;; register.
  2285. ;;;
  2286. #+32x16-divide
  2287. (defun 32x16-multiply (a b)
  2288.   (declare (type (unsigned-byte 16) a b))
  2289.   (the bignum-element-type (* a b)))
  2290.  
  2291. ;;; 32X16-MULTIPLY-SPLIT -- Internal.        [optionally IN ASSEMBLER]
  2292. ;;;
  2293. ;;; This multiplies a and b, 16-bit quantities, and returns the result as two
  2294. ;;; 16-bit quantities, high and low.
  2295. ;;;
  2296. #+32x16-divide
  2297. (defun 32x16-multiply-split (a b)
  2298.   (let ((res (32x16-multiply a b)))
  2299.     (declare (the bignum-element-type res))
  2300.     (values (the (unsigned-byte 16) (logand #xFFFF (ash res -16)))
  2301.         (the (unsigned-byte 16) (logand #xFFFF res)))))
  2302.  
  2303.  
  2304.  
  2305. ;;; The %FLOOR below uses this buffer the same way BIGNUM-TRUNCATE uses
  2306. ;;; *truncate-x*.  There's no y buffer since we pass around the two 16-bit
  2307. ;;; digits and use them slightly differently than the general truncation
  2308. ;;; algorithm above.
  2309. ;;;
  2310. #+32x16-divide
  2311. (defvar *32x16-truncate-x* (make-array 4 :element-type '(unsigned-byte 16)
  2312.                        :initial-element 0))
  2313.  
  2314. ;;; %FLOOR -- Internal.        LEFT IMPLEMENTED AT LISP LEVEL
  2315. ;;;
  2316. ;;; This does the same thing as the %FLOOR above, but it does it at Lisp level
  2317. ;;; when there is no 64x32-bit divide instruction on the machine.
  2318. ;;;
  2319. ;;; It implements the higher level tactics of BIGNUM-TRUNCATE, but it makes use
  2320. ;;; of special situation provided, four 16-bit digits divided by two 16-bit
  2321. ;;; digits.
  2322. ;;;
  2323. #+32x16-divide
  2324. (defun %floor (a b c)
  2325.   (declare (type bignum-element-type a b c))
  2326.   ;;
  2327.   ;; Setup *32x16-truncate-x* buffer from a and b.
  2328.   (setf (aref *32x16-truncate-x* 0)
  2329.     (the (unsigned-byte 16) (logand #xFFFF b)))
  2330.   (setf (aref *32x16-truncate-x* 1)
  2331.     (the (unsigned-byte 16)
  2332.          (logand #xFFFF
  2333.              (the (unsigned-byte 16) (ash b -16)))))
  2334.   (setf (aref *32x16-truncate-x* 2)
  2335.     (the (unsigned-byte 16) (logand #xFFFF a)))
  2336.   (setf (aref *32x16-truncate-x* 3)
  2337.     (the (unsigned-byte 16)
  2338.          (logand #xFFFF
  2339.              (the (unsigned-byte 16) (ash a -16)))))
  2340.   ;;
  2341.   ;; From DO-TRUNCATE, but unroll the loop.
  2342.   (let* ((y1 (logand #xFFFF (ash c -16)))
  2343.      (y2 (logand #xFFFF c))
  2344.      (q (the bignum-element-type
  2345.          (ash (32x16-try-bignum-truncate-guess
  2346.                (32x16-truncate-guess y1 y2
  2347.                          (aref *32x16-truncate-x* 3)
  2348.                          (aref *32x16-truncate-x* 2)
  2349.                          (aref *32x16-truncate-x* 1))
  2350.                y1 y2 1)
  2351.               16))))
  2352.     (declare (type bignum-element-type q)
  2353.          (type (unsigned-byte 16) y1 y2))
  2354.     (values (the bignum-element-type
  2355.          (logior q
  2356.              (the (unsigned-byte 16)
  2357.                   (32x16-try-bignum-truncate-guess
  2358.                    (32x16-truncate-guess
  2359.                 y1 y2
  2360.                 (aref *32x16-truncate-x* 2)
  2361.                 (aref *32x16-truncate-x* 1)
  2362.                 (aref *32x16-truncate-x* 0))
  2363.                    y1 y2 0))))
  2364.         (the bignum-element-type
  2365.          (logior (the bignum-element-type
  2366.                   (ash (aref *32x16-truncate-x* 1) 16))
  2367.              (the (unsigned-byte 16)
  2368.                   (aref *32x16-truncate-x* 0)))))))
  2369.  
  2370. ;;; 32X16-TRY-BIGNUM-TRUNCATE-GUESS  --  Internal.
  2371. ;;;
  2372. ;;; This is similar to TRY-BIGNUM-TRUNCATE-GUESS, but this unrolls the two
  2373. ;;; loops.  This also substitutes for %DIGIT-0-OR-PLUSP the equivalent
  2374. ;;; expression without any embellishment or pretense of abstraction.  The first
  2375. ;;; loop is unrolled, but we've put the body of the loop into the function
  2376. ;;; 32X16-TRY-GUESS-ONE-RESULT-DIGIT.
  2377. ;;;
  2378. #+32x16-divide
  2379. (defun 32x16-try-bignum-truncate-guess (guess y-high y-low low-x-digit)
  2380.   (declare (type bignum-index low-x-digit)
  2381.        (type (unsigned-byte 16) guess y-high y-low))
  2382.   (let ((high-x-digit (+ 2 low-x-digit)))
  2383.     ;;
  2384.     ;; Multiply guess and divisor, subtracting from dividend simultaneously.
  2385.     (multiple-value-bind
  2386.     (guess*y-hold carry borrow)
  2387.     (32x16-try-guess-one-result-digit guess y-low 0 0 1 low-x-digit)
  2388.       (declare (type (unsigned-byte 16) guess*y-hold)
  2389.            (fixnum carry borrow))
  2390.       (multiple-value-bind
  2391.       (guess*y-hold carry borrow)
  2392.       (32x16-try-guess-one-result-digit guess y-high guess*y-hold
  2393.                         carry borrow (1+ low-x-digit))
  2394.     (declare (type (unsigned-byte 16) guess*y-hold)
  2395.          (fixnum borrow)
  2396.          (ignore carry))
  2397.     (setf (aref *32x16-truncate-x* high-x-digit)
  2398.           (32x16-subtract-with-borrow (aref *32x16-truncate-x* high-x-digit)
  2399.                       guess*y-hold borrow))))
  2400.     ;;
  2401.     ;; See if guess is off by one, adding one Y back in if necessary.
  2402.     (cond ((zerop (logand #x8000 (aref *32x16-truncate-x* high-x-digit)))
  2403.        ;; The subtraction result is zero or positive.
  2404.        guess)
  2405.       (t
  2406.        ;; If subtraction has negative result, add one divisor value back in.
  2407.        ;; The guess was one two large in magnitude.
  2408.        (multiple-value-bind (v carry)
  2409.                 (32x16-add-with-carry y-low
  2410.                               (aref *32x16-truncate-x*
  2411.                                 low-x-digit)
  2412.                               0)
  2413.          (declare (type (unsigned-byte 16) v))
  2414.          (setf (aref *32x16-truncate-x* low-x-digit) v)
  2415.          (multiple-value-bind (v carry)
  2416.                   (32x16-add-with-carry y-high
  2417.                             (aref *32x16-truncate-x*
  2418.                                   (1+ low-x-digit))
  2419.                             carry)
  2420.            (setf (aref *32x16-truncate-x* (1+ low-x-digit)) v)
  2421.            (setf (aref *32x16-truncate-x* high-x-digit)
  2422.              (32x16-add-with-carry (aref *32x16-truncate-x* high-x-digit)
  2423.                        carry 0))))
  2424.        (if (zerop (logand #x8000 guess))
  2425.            (1- guess)
  2426.            (1+ guess))))))
  2427.  
  2428. ;;; 32X16-TRY-GUESS-ONE-RESULT-DIGIT -- Internal.
  2429. ;;;
  2430. ;;; This is similar to the body of the loop in TRY-BIGNUM-TRUNCATE-GUESS that
  2431. ;;; multiplies the guess by y and subtracts the result from x simultaneously.
  2432. ;;; This returns the digit remembered as part of the multiplication, the carry
  2433. ;;; from additions done on behalf of the multiplication, and the borrow from
  2434. ;;; doing the subtraction.
  2435. ;;;
  2436. #+32x16-divide
  2437. (defun 32x16-try-guess-one-result-digit (guess y-digit guess*y-hold
  2438.                      carry borrow x-index)
  2439.   (multiple-value-bind (high-digit low-digit)
  2440.                (32x16-multiply-split guess y-digit)
  2441.     (declare (type (unsigned-byte 16) high-digit low-digit))
  2442.     (multiple-value-bind (low-digit temp-carry)
  2443.              (32x16-add-with-carry low-digit guess*y-hold carry)
  2444.       (declare (type (unsigned-byte 16) low-digit))
  2445.       (multiple-value-bind (high-digit temp-carry)
  2446.                (32x16-add-with-carry high-digit temp-carry 0)
  2447.     (declare (type (unsigned-byte 16) high-digit))
  2448.     (multiple-value-bind (x temp-borrow)
  2449.                  (32x16-subtract-with-borrow
  2450.                   (aref *32x16-truncate-x* x-index)
  2451.                   low-digit borrow)
  2452.       (declare (type (unsigned-byte 16) x))
  2453.       (setf (aref *32x16-truncate-x* x-index) x)
  2454.       (values high-digit temp-carry temp-borrow))))))
  2455.  
  2456. ;;; 32X16-TRUNCATE-GUESS -- Internal.
  2457. ;;;
  2458. ;;; This is similar to BIGNUM-TRUNCATE-GUESS, but instead of computing the
  2459. ;;; guess exactly as described in the its comments (digit by digit), this
  2460. ;;; massages the 16-bit quantities into 32-bit quantities and performs the
  2461. ;;; 
  2462. #+32x16-divide
  2463. (defun 32x16-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
  2464.   (declare (type (unsigned-byte 16) y1 y2 x-i x-i-1 x-i-2))
  2465.   (let ((guess (if (= x-i y1)
  2466.            #xFFFF
  2467.            (32x16-divide x-i x-i-1 y1))))
  2468.     (declare (type (unsigned-byte 16) guess))
  2469.     (loop
  2470.       (let* ((guess*y1 (the bignum-element-type
  2471.                 (ash (logand #xFFFF
  2472.                      (the bignum-element-type
  2473.                           (32x16-multiply guess y1)))
  2474.                  16)))
  2475.          (x-y (%subtract-with-borrow
  2476.            (the bignum-element-type
  2477.             (logior (the bignum-element-type
  2478.                      (ash x-i-1 16))
  2479.                 x-i-2))
  2480.            guess*y1
  2481.            1))
  2482.          (guess*y2 (the bignum-element-type (%multiply guess y2))))
  2483.     (declare (type bignum-element-type guess*y1 x-y guess*y2))
  2484.     (if (%digit-greater guess*y2 x-y)
  2485.         (decf guess)
  2486.         (return guess))))))
  2487.  
  2488.  
  2489. ;;;; General utilities.
  2490.  
  2491. ;;; MAKE-SMALL-BIGNUM -- Public.
  2492. ;;;
  2493. ;;; Allocate a single word bignum that holds fixnum.  This is useful when
  2494. ;;; we are trying to mix fixnum and bignum operands.
  2495. ;;; 
  2496. (proclaim '(inline make-small-bignum))
  2497. (defun make-small-bignum (fixnum)
  2498.   (let ((res (%allocate-bignum 1)))
  2499.     (setf (%bignum-ref res 0) (%fixnum-to-digit fixnum))
  2500.     res))
  2501.  
  2502. ;;; %NORMALIZE-BIGNUM-BUFFER -- Internal.
  2503. ;;;
  2504. ;;; Internal in-place operations use this to fixup remaining digits in the
  2505. ;;; incoming data, such as in-place shifting.  This is basically the same as
  2506. ;;; the first form in %NORMALIZE-BIGNUM, but we return the length of the buffer
  2507. ;;; instead of shrinking the bignum.
  2508. ;;;
  2509. #+nil(proclaim '(ext:maybe-inline %normalize-bignum-buffer))
  2510. (defun %normalize-bignum-buffer (result len)
  2511.   (declare (type bignum-type result)
  2512.        (type bignum-index len))
  2513.   (unless (= len 1)
  2514.     (do ((next-digit (%bignum-ref result (- len 2))
  2515.              (%bignum-ref result (- len 2)))
  2516.      (sign-digit (%bignum-ref result (1- len)) next-digit))
  2517.     ((not (zerop (logxor sign-digit (%ashr next-digit (1- digit-size))))))
  2518.       (when (= (decf len) 1)
  2519.     (return))
  2520.       (setf (%bignum-ref result len) 0)))
  2521.   len)
  2522.  
  2523. ;;; %NORMALIZE-BIGNUM -- Internal.
  2524. ;;;
  2525. ;;; This drops the last digit if it is unnecessary sign information.  It
  2526. ;;; repeats this as needed, possibly ending with a fixnum.  If the resulting
  2527. ;;; length from shrinking is one, see if our one word is a fixnum.  Shift the
  2528. ;;; possible fixnum bits completely out of the word, and compare this with
  2529. ;;; shifting the sign bit all the way through.  If the bits are all 1's or 0's
  2530. ;;; in both words, then there are just sign bits between the fixnum bits and
  2531. ;;; the sign bit.  If we do have a fixnum, shift it over for the two low-tag
  2532. ;;; bits.
  2533. ;;;
  2534. (defun %normalize-bignum (result len)
  2535.   (declare (type bignum-type result)
  2536.        (type bignum-index len)
  2537.        #+nil(inline %normalize-bignum-buffer))
  2538.   (let ((newlen (%normalize-bignum-buffer result len)))
  2539.     (declare (type bignum-index newlen))
  2540.     (unless (= newlen len)
  2541.       (%bignum-set-length result newlen))
  2542.     (if (= newlen 1)
  2543.     (let ((digit (%bignum-ref result 0)))
  2544.       (if (= (%ashr digit 29) (%ashr digit (1- digit-size)))
  2545.           (%fixnum-digit-with-correct-sign digit)
  2546.           result))
  2547.     result)))
  2548.  
  2549. ;;; %MOSTLY-NORMALIZE-BIGNUM -- Internal.
  2550. ;;;
  2551. ;;; This drops the last digit if it is unnecessary sign information.  It
  2552. ;;; repeats this as needed, possibly ending with a fixnum magnitude but never
  2553. ;;; returning a fixnum.
  2554. ;;;
  2555. (defun %mostly-normalize-bignum (result len)
  2556.   (declare (type bignum-type result)
  2557.        (type bignum-index len)
  2558.        #+nil(inline %normalize-bignum-buffer))
  2559.   (let ((newlen (%normalize-bignum-buffer result len)))
  2560.     (declare (type bignum-index newlen))
  2561.     (unless (= newlen len)
  2562.       (%bignum-set-length result newlen))
  2563.     result))
  2564.