Skip to content
Snippets Groups Projects
bignum.lisp 127 KiB
Newer Older
;;; -*- Mode: Lisp; Log: code.log; Package: bignum -*-
wlott's avatar
wlott committed
;;;
;;; **********************************************************************
;;; This code was written as part of the CMU Common Lisp project at
;;; Carnegie Mellon University, and has been placed in the public domain.
;;;
(ext:file-comment
  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/bignum.lisp,v 1.40 2004/08/31 00:19:41 rtoy Exp $")
wlott's avatar
wlott committed
;;; **********************************************************************
;;;
;;; This file contains code to implement bignum support.
;;;
wlott's avatar
wlott committed
(in-package "BIGNUM")
(use-package "KERNEL")

;;; These symbols define the interface to the number code.

(export '(add-bignums multiply-bignums negate-bignum subtract-bignum
	  multiply-bignum-and-fixnum multiply-fixnums
	  bignum-ashift-right bignum-ashift-left bignum-gcd
	  bignum-to-float bignum-integer-length
	  bignum-logical-and bignum-logical-ior bignum-logical-xor
	  bignum-logical-not bignum-load-byte bignum-deposit-byte
	  bignum-truncate bignum-plus-p bignum-compare make-small-bignum
	  bignum-logcount bignum-logbitp))

;;; These symbols define the interface to the compiler.

(export '(bignum-type bignum-element-type bignum-index %allocate-bignum
	  %bignum-length %bignum-set-length %bignum-ref %bignum-set
	  %digit-0-or-plusp %add-with-carry %subtract-with-borrow
	  %multiply-and-add %multiply %lognot %logand %logior %logxor
	  %fixnum-to-digit %floor %fixnum-digit-with-correct-sign %ashl
	  %ashr %digit-logical-shift-right))

wlott's avatar
wlott committed


;;;; Notes.

;;; The following interfaces will either be assembler routines or code sequences
;;; expanded into the code as basic bignum operations:
;;;    General:
;;;       %BIGNUM-LENGTH
;;;       %ALLOCATE-BIGNUM
;;;       %BIGNUM-REF
;;;       %NORMALIZE-BIGNUM
wlott's avatar
wlott committed
;;;       %FIXNUM-DIGIT-WITH-CORRECT-SIGN
;;;       %SIGN-DIGIT
;;;	  %ASHR
;;;       %ASHL
;;;       %BIGNUM-0-OR-PLUSP
;;;       %DIGIT-LOGICAL-SHIFT-RIGHT
wlott's avatar
wlott committed
;;;    General (May not exist when done due to sole use in %-routines.)
;;;       %DIGIT-0-OR-PLUSP
;;;    Addition:
;;;       %ADD-WITH-CARRY
;;;    Subtraction:
;;;       %SUBTRACT-WITH-BORROW
;;;    Multiplication
;;;       %MULTIPLY
;;;    Negation
;;;       %LOGNOT
;;;    Shifting (in place)
;;;       %NORMALIZE-BIGNUM-BUFFER
;;;    GCD/Relational operators:
;;;       %DIGIT-COMPARE
;;;       %DIGIT-GREATER
wlott's avatar
wlott committed
;;;    Relational operators:
;;;       %LOGAND
;;;       %LOGIOR
;;;       %LOGXOR
;;;    LDB
;;;       %FIXNUM-TO-DIGIT
;;;    TRUNCATE
;;;       %FLOOR
;;;
;;;
;;; Note: The floating routines know about the float representation.
;;;
wlott's avatar
wlott committed
;;; PROBLEM 1:
;;; There might be a problem with various LET's and parameters that take a
;;; digit value.  We need to write these so those things stay in 32-bit
;;; registers and number stack slots.  I bind locals to these values, and I
;;; use function on them -- ZEROP, ASH, etc.
;;;
;;; PROBLEM 2:
;;; In shifting and byte operations, I use masks and logical operations that
;;; could result in intermediate bignums.  This is hidden by the current system,
;;; but I may need to write these in a way that keeps these masks and logical
;;; operations from diving into the Lisp level bignum code.
;;;
;;; To do:
;;;    fixnums
;;;       logior, logxor, logand
;;;       depending on relationals, < (twice) and <= (twice)
;;;          or write compare thing (twice).
;;;       LDB on fixnum with bignum result.
;;;       DPB on fixnum with bignum result.
;;;       TRUNCATE returns zero or one as one value and fixnum or minus fixnum
;;;          for the other value when given (truncate fixnum bignum).
;;;          Returns (truncate bignum fixnum) otherwise.
;;;       addition
;;;       subtraction (twice)
;;;       multiply
;;;       GCD
;;;    write MASK-FIELD and DEPOSIT-FIELD in terms of logical operations.
;;;    DIVIDE
;;;       IF (/ x y) with bignums:
;;;          do the truncate, and if rem is 0, return quotient.
;;;          if rem is non-0
;;;	     gcd of x and y.
;;;	     "truncate" each by gcd, ignoring remainder 0.
;;;	     form ratio of each result, bottom is positive.
;;;



;;;; What's a bignum?

(eval-when (compile load eval) ;Necessary for DEFTYPE.

cwang's avatar
cwang committed
(defconstant digit-size #+amd64 64 #-amd64 32)
wlott's avatar
wlott committed

(defconstant maximum-bignum-length (1- (ash 1 (- 32 vm:type-bits))))
wlott's avatar
wlott committed

) ;eval-when



;;;; Internal inline routines.

;;; %ALLOCATE-BIGNUM must zero all elements.
;;;
(defun %allocate-bignum (length)
  (declare (type bignum-index length))
  (%allocate-bignum length))

;;; Extract the length of the bignum.
;;; 
(defun %bignum-length (bignum)
  (declare (type bignum-type bignum))
  (%bignum-length bignum))

;;; %BIGNUM-REF needs to access bignums as obviously as possible, and it needs
;;; to be able to return 32 bits somewhere no one looks for real objects.
;;;
(defun %bignum-ref (bignum i)
  (declare (type bignum-type bignum)
	   (type bignum-index i))
  (%bignum-ref bignum i))
;;;
(defun %bignum-set (bignum i value)
  (declare (type bignum-type bignum)
	   (type bignum-index i)
	   (type bignum-element-type value))
  (%bignum-set bignum i value))
;;;
(defsetf %bignum-ref %bignum-set)

;;; Return T if digit is positive, or NIL if negative.
;;;
(defun %digit-0-or-plusp (digit)
  (declare (type bignum-element-type digit))
  (not (logbitp (1- digit-size) digit)))
wlott's avatar
wlott committed

pw's avatar
pw committed
(declaim (inline %bignum-0-or-plusp))
wlott's avatar
wlott committed
(defun %bignum-0-or-plusp (bignum len)
  (declare (type bignum-type bignum)
	   (type bignum-index len))
  (%digit-0-or-plusp (%bignum-ref bignum (1- len))))

;;; %ADD-WITH-CARRY -- Internal.
;;;
;;; This should be in assembler, and should not cons intermediate results.  It
;;; returns a 32bit digit and a carry resulting from adding together a, b, and
;;; an incoming carry.
;;;
(defun %add-with-carry (a b carry)
  (declare (type bignum-element-type a b)
	   (type (mod 2) carry))
  (%add-with-carry a b carry))

;;; %SUBTRACT-WITH-BORROW -- Internal.
;;;
;;; This should be in assembler, and should not cons intermediate results.  It
;;; returns a 32bit digit and a borrow resulting from subtracting b from a, and
;;; subtracting a possible incoming borrow.
;;;
;;; We really do:  a - b - 1 + borrow, where borrow is either 0 or 1.
;;; 
(defun %subtract-with-borrow (a b borrow)
  (declare (type bignum-element-type a b)
	   (type (mod 2) borrow))
  (%subtract-with-borrow a b borrow))

;;; %MULTIPLY -- Internal.
;;;
;;; This multiplies two digit-size (32-bit) numbers, returning a 64-bit result
;;; split into two 32-bit quantities.
;;;
(defun %multiply (x y)
  (declare (type bignum-element-type x y))
  (%multiply x y))

;;; %MULTIPLY-AND-ADD  --  Internal.
;;;
;;; This multiplies x-digit and y-digit, producing high and low digits
;;; manifesting the result.  Then it adds the low digit, res-digit, and
;;; carry-in-digit.  Any carries (note, you still have to add two digits at a
;;; time possibly producing two carries) from adding these three digits get
;;; added to the high digit from the multiply, producing the next carry digit.
;;; Res-digit is optional since two uses of this primitive multiplies a single
;;; digit bignum by a multiple digit bignum, and in this situation there is no
;;; need for a result buffer accumulating partial results which is where the
;;; res-digit comes from.
;;;
(defun %multiply-and-add (x-digit y-digit carry-in-digit &optional (res-digit 0))
  (declare (type bignum-element-type x-digit y-digit res-digit carry-in-digit))
  (%multiply-and-add x-digit y-digit carry-in-digit res-digit))
wlott's avatar
wlott committed
;;; %LOGNOT -- Internal.
;;;
(defun %lognot (digit)
  (declare (type bignum-element-type digit))
  (%lognot digit))

;;; %LOGAND -- Internal.
;;; %LOGIOR -- Internal.
;;; %LOGXOR -- Internal.
;;;
;;; Do the 32bit unsigned op.
;;;
pw's avatar
pw committed
(declaim (inline %logand %logior %logxor))
wlott's avatar
wlott committed
(defun %logand (a b)
  (declare (type bignum-element-type a b))
  (logand a b))
(defun %logior (a b)
  (declare (type bignum-element-type a b))
  (logior a b))
(defun %logxor (a b)
  (declare (type bignum-element-type a b))
  (logxor a b))

;;; %FIXNUM-TO-DIGIT -- Internal.
;;;
;;; This takes a fixnum and sets it up as an unsigned 32-bit quantity.  In
;;; the new system this will mean shifting it right two bits.
;;;
(defun %fixnum-to-digit (x)
  (declare (fixnum x))
  (logand x (1- (ash 1 digit-size))))

wlott's avatar
wlott committed
;;; %FLOOR -- Internal.
;;;
;;; This takes three digits and returns the FLOOR'ed result of dividing the
;;; first two as a 64-bit integer by the third.
;;;
;;; DO WEIRD let AND setq STUFF TO SLIME THE COMPILER INTO ALLOWING THE %FLOOR
;;; TRANSFORM TO EXPAND INTO PSEUDO-ASSEMBLER FOR WHICH THE COMPILER CAN LATER
;;; CORRECTLY ALLOCATE REGISTERS.
;;;
wlott's avatar
wlott committed
(defun %floor (a b c)
  (let ((a a) (b b) (c c))
    (declare (type bignum-element-type a b c))
wlott's avatar
wlott committed


;;; %FIXNUM-DIGIT-WITH-CORRECT-SIGN -- Internal.
;;;
;;; Convert the digit to a regular integer assuming that the digit is signed.
;;;
(defun %fixnum-digit-with-correct-sign (digit)
  (declare (type bignum-element-type digit))
  (if (logbitp (1- digit-size) digit)
      (logior digit (ash -1 digit-size))
      digit))

;;; %ASHR -- Internal.
;;;
;;; Do an arithmetic shift right of data even though bignum-element-type is
;;; unsigned.
;;;
(defun %ashr (data count)
  (declare (type bignum-element-type data)
cwang's avatar
cwang committed
	   (type (mod #+amd64 64 #-amd64 32) count))
wlott's avatar
wlott committed
  (%ashr data count))

;;; %ASHL -- Internal.
;;;
;;; This takes a 32-bit quantity and shifts it to the left, returning a 32-bit
;;; quantity.
(defun %ashl (data count)
  (declare (type bignum-element-type data)
	   (type (mod 32) count))
  (%ashl data count))

;;; %DIGIT-LOGICAL-SHIFT-RIGHT  --  Internal
;;;
;;;    Do an unsigned (logical) right shift of a digit by Count.
;;;
(defun %digit-logical-shift-right (data count)
  (declare (type bignum-element-type data)
	   (type (mod 32) count))
  (%digit-logical-shift-right data count))


wlott's avatar
wlott committed
;;; %BIGNUM-SET-LENGTH -- Internal.
;;;
;;; Change the length of bignum to be newlen.  Newlen must be the same or
;;; smaller than the old length, and any elements beyond newlen must be zeroed.
;;; 
(defun %bignum-set-length (bignum newlen)
  (declare (type bignum-type bignum)
	   (type bignum-index newlen))
  (%bignum-set-length bignum newlen))

;;; %SIGN-DIGIT -- Internal.
;;;
;;; This returns 0 or "-1" depending on whether the bignum is positive.  This
;;; is suitable for infinite sign extension to complete additions,
;;; subtractions, negations, etc.  This cannot return a -1 represented as
;;; a negative fixnum since it would then have to low zeros.
;;;
pw's avatar
pw committed
(declaim (inline %sign-digit))
wlott's avatar
wlott committed
(defun %sign-digit (bignum len)
  (declare (type bignum-type bignum)
	   (type bignum-index len))
  (%ashr (%bignum-ref bignum (1- len)) (1- digit-size)))



;;; %DIGIT-COMPARE and %DIGIT-GREATER -- Internal.
;;;
;;; These take two 32 bit quantities and compare or contrast them without
;;; wasting time with incorrect type checking.
;;;
pw's avatar
pw committed
(declaim (inline %digit-compare %digit-greater))
(defun %digit-compare (x y)
  (= x y))
;;;
(defun %digit-greater (x y)
  (> x y))

wlott's avatar
wlott committed

(declaim (optimize (speed 3) (safety 0)(ext:inhibit-warnings 3)))
wlott's avatar
wlott committed


;;;; Addition.

(defun add-bignums (a b)
  (declare (type bignum-type a b))
  (let ((len-a (%bignum-length a))
	(len-b (%bignum-length b)))
    (declare (type bignum-index len-a len-b))
    (multiple-value-bind (a len-a b len-b)
			 (if (> len-a len-b)
			     (values a len-a b len-b)
			     (values b len-b a len-a))
      (declare (type bignum-type a b)
	       (type bignum-index len-a len-b))
      (let* ((len-res (1+ len-a))
	     (res (%allocate-bignum len-res))
	     (carry 0))
	(declare (type bignum-index len-res)
		 (type bignum-type res)
		 (type (mod 2) carry))
	(dotimes (i len-b)
	  (declare (type bignum-index i))
	  (multiple-value-bind
	      (v k)
	      (%add-with-carry (%bignum-ref a i) (%bignum-ref b i) carry)
	    (declare (type bignum-element-type v)
		     (type (mod 2) k))
	    (setf (%bignum-ref res i) v)
	    (setf carry k)))
	(if (/= len-a len-b)
	    (finish-add a res carry (%sign-digit b len-b) len-b len-a)
	    (setf (%bignum-ref res len-a)
		  (%add-with-carry (%sign-digit a len-a)
				   (%sign-digit b len-b)
				   carry)))
	(%normalize-bignum res len-res)))))

;;; FINISH-ADD -- Internal.
;;;
;;; This takes the longer of two bignums and propagates the carry through its
;;; remaining high order digits.
;;;
(defun finish-add (a res carry sign-digit-b start end)
  (declare (type bignum-type a res)
	   (type (mod 2) carry)
	   (type bignum-element-type sign-digit-b)
	   (type bignum-index start end))
  (do ((i start (1+ i)))
      ((= i end)
       (setf (%bignum-ref res end)
	     (%add-with-carry (%sign-digit a end) sign-digit-b carry)))
    (declare (type bignum-index i))
wlott's avatar
wlott committed
    (multiple-value-bind (v k)
			 (%add-with-carry (%bignum-ref a i) sign-digit-b carry)
      (setf (%bignum-ref res i) v)
      (setf carry k)))
  (ext:undefined-value))
wlott's avatar
wlott committed


;;;; Subtraction.

(eval-when (compile eval)

;;; SUBTRACT-BIGNUM-LOOP -- Internal.
;;;
;;; This subtracts b from a plugging result into res.  Return-fun is the
;;; function to call that fixes up the result returning any useful values, such
;;; as the result.  This macro may evaluate its arguments more than once.
;;;
(defmacro subtract-bignum-loop (a len-a b len-b res len-res return-fun)
  (let ((borrow (gensym))
	(a-digit (gensym))
	(a-sign (gensym))
	(b-digit (gensym))
	(b-sign (gensym))
wlott's avatar
wlott committed
	(i (gensym))
	(v (gensym))
	(k (gensym)))
    `(let* ((,borrow 1)
	    (,a-sign (%sign-digit ,a ,len-a))
	    (,b-sign (%sign-digit ,b ,len-b)))
       (declare (type bignum-element-type ,a-sign ,b-sign))
       (dotimes (,i ,len-res)
	 (declare (type bignum-index ,i))
	 (let ((,a-digit (if (< ,i ,len-a) (%bignum-ref ,a ,i) ,a-sign))
	       (,b-digit (if (< ,i ,len-b) (%bignum-ref ,b ,i) ,b-sign)))
	   (declare (type bignum-element-type ,a-digit ,b-digit))
	   (multiple-value-bind
	       (,v ,k)
	       (%subtract-with-borrow ,a-digit ,b-digit ,borrow)
	     (setf (%bignum-ref ,res ,i) ,v)
	     (setf ,borrow ,k))))
wlott's avatar
wlott committed
       (,return-fun ,res ,len-res))))

) ;EVAL-WHEN

(defun subtract-bignum (a b)
  (declare (type bignum-type a b))
  (let* ((len-a (%bignum-length a))
	 (len-b (%bignum-length b))
	 (len-res (1+ (max len-a len-b)))
wlott's avatar
wlott committed
	 (res (%allocate-bignum len-res)))
    (declare (type bignum-index len-a len-b len-res)) ;Test len-res for bounds?
    (subtract-bignum-loop a len-a b len-b res len-res %normalize-bignum)))

;;; SUBTRACT-BIGNUM-BUFFERS -- Internal.
;;;
;;; Operations requiring a subtraction without the overhead of intermediate
;;; results, such as GCD, use this.  It assumes Result is big enough for the
;;; result.
;;;
(defun subtract-bignum-buffers-with-len (a len-a b len-b result len-res)
wlott's avatar
wlott committed
  (declare (type bignum-type a b)
	   (type bignum-index len-a len-b))
  (subtract-bignum-loop a len-a b len-b result len-res
			%normalize-bignum-buffer))

(defun subtract-bignum-buffers (a len-a b len-b result)
  (declare (type bignum-type a b)
	   (type bignum-index len-a len-b))  
  (subtract-bignum-loop a len-a b len-b result (max len-a len-b)
			%normalize-bignum-buffer))
wlott's avatar
wlott committed



;;;; Multiplication.

;;; **********************************************************************
;;;
;;; This Karatsuba implementation was written by Raymond Toy and
;;; Douglas T. Crosher and has been placed in the Public domain, and
;;; is provided 'as is'.
;;;
;;; **********************************************************************
;;;
;;; Karatsuba multiplication.
;;; References:  Knuth, Seminumerical Algorithms.
;;;

;;; Based on an implementation by Raymond Toy, who's explanation follows:
;;;
;;; Let U = (b^n)*U1 + U0 and V = (b^n)*V1 + V0 be two positive
;;; integers consisting of b^(2n) base-b digits.  Typically, we use
;;; either base 2 digits or base 2^32 digits.
;;;
;;; The product U*V can be written
;;;
;;;    U*V = (b^(2n))*U1*V1 + (b^n)(U0*V1 + U1*V0) + U0*V0
;;;
;;; However, the middle term can be written as
;;;
;;;    U0*V1 + U1*V0 = U0*V0 + U1*V1 - (U0 - U1)*(V0 - V1) 
;;;
;;; Thus,
;;;
;;;    U*V = (b^(2n))*U1+V1 + (b^n)*[U1*V1 + U0*V0 - (U0 - U1)*(V0 - V1)]
;;;           U0*V0
;;;
;;; Note that now we require only three multiplies: U1*V1, U0*V0, and
;;; (U0 - U1)*(V0 - V1).  This should speed up multiplication
;;; considerably.
;;;
;;; For implemenation, we only want to work with positive numbers.
;;; Since U0-U1 or V0-V1 may be negative, we need to account for this.
;;; To do so, we compute
;;;
;;;    U1*V1 + U0*V0 - Su*Sv*|U0 - U1|*|V0 - V1|
;;;
;;; where Su = sign (U0 - U1) and Sv = sign (V0 - V1).
;;;
;;; Note that some authors use the expression
;;;
;;;    U0*V1 + U1*V0 = (U0 + U1)*(V0 + V1) - U0*V0 - U1*V1
;;;
;;; We don't because if U0 and U1 each contain b^n digits, U0 + U1 can
;;; overflow to b^n + 1 digits.  This complicates the multiplication quite a
;;; bit, so we use |U0 - U1|.  This is guaranteed not to overflow since U0 and
;;; U1 are both positive with b^n digits.
;;;

;;; karatsuba-internal  --  Internal.
;;;
;;; Main routine for Karatsuba multiplication.
;;; By Douglas Crosher.
;;;
;;; The result is placed in the result bignum starting at position RET, and is
;;; of size 2xD words. An extra D words are used during the calculations.
;;;
;;; The calculations are packed into this working space as follows:
;;;
;;; ------------------------------------------------------------------
;;; |    |   D/2   |   D/2   |  D/2    |   D/2   |   D/2   |  D/2    |
;;; |----+---------+---------+---------+---------+---------+---------|
;;; | A. | |u1-u0| | |v1-v0| |    x    |       |u1-u0|*|v1-v0|       | 
;;; | B. |           u0 * v0           |       |u1-u0|*|v1-v0|       |
;;; | C. | u0 * v0 - S|u1-u0|*|v1-v0|  |             x               |
;;; | D. | u0 * v0 - S|u1-u0|*|v1-v0|  |          v1 * v0            |
;;; | E. | u0 * v0 - S|u1-u0|*|v1-v0| + v1 * v0  |    x    |    x    |
;;; ------------------------------------------------------------------
;;;
;;; A. Calculate |u1-u0| place the result in the first D/2 words, and
;;;    calculate |v1-v0| placing the result the second D/2 words. Then
;;;    calculate |u1-u0|*|v1-v0| using the last 3 D/2 words for the result and
;;;    subcalculate scratch area.
;;;
;;; B. Calculate the product u0*v0, using the first 3 D/2 words for the result
;;;    and subcalculate scratch area.
;;;
;;; C. Sum the results of first two stages into the first 3 D/2 words.
;;;    The calculation is (u0*v0) + ((u0*v0) - S|u1-u0|*|v1-v0|) *2^(d/2).
;;;
;;; D. Calculate the product u1*v1, using the last 3 D/2 words for the result
;;;    and subcalculate scratch area.
;;;
;;; E. Finally sum in u1*v1 to give the result in the first 4 D/2 words.
;;;


;; This might need some tuning
(declaim (fixnum *karatsuba-classical-cutoff*))
(defparameter *karatsuba-classical-cutoff* 10
  "When the bignum pieces are smaller than this many words, we use the
classical multiplication algorithm instead of recursing all the way
down to individual words.")

(defun karatsuba-internal (u v ret n d result)
  (declare (type bignum-index u v ret)
	   (type bignum-index d n)
	   (type bignum-type result)
	   (optimize (speed 3) (safety 0) (debug 0)))

  (let* ((odd (logand d 1))
	 (origd d)
	 (d/2- (ash d -1))
	 (d/2 (+ d/2- odd))
	 (d (* d/2 2))
	 (n/2 (ash n -1)))
    (labels ((bignum-multiply-in-place (u v ret d result cleard)
	       (declare (type bignum-index u v ret)
			(type bignum-index d)
			(type bignum-type result)
			(optimize (speed 3) (safety 0) (debug 0)))
	       ;; Clear out the result area, which might contain garbage for
	       ;; other recursions!
	       (do ((i ret (1+ i))
		    (end (+ ret (* 2 cleard))))
		   ((>= i end))
		 (setf (%bignum-ref result i) 0))
	       ;; Stolen from bignum.lisp multiply-bignum: a straight forward
	       ;; multiplication via the classical schoolbook algorithm.
	       (dotimes (i d)
		 (declare (type bignum-index i))
		 (let ((carry-digit 0)
		       (x (%bignum-ref result (+ u i)))
		       (k (+ ret i)))
		   (declare (type bignum-index k)
			    (type bignum-element-type carry-digit x))
		   (do ((j d (1- j))
			(v+j v (1+ v+j)))
		       ((zerop j))
		     (declare (type bignum-index j v+j))
		     (multiple-value-bind (big-carry res-digit)
			 (%multiply-and-add x
					    (%bignum-ref result v+j)
					    (%bignum-ref result k)
					    carry-digit)
		       (declare (type bignum-element-type big-carry res-digit))
		       (setf (%bignum-ref result k) res-digit)
		       (setf carry-digit big-carry)
		       (incf k)))
		   (setf (%bignum-ref result k) carry-digit))))
	     (abs-diff (n1 x n2 y res result)
	       (declare (type bignum-index n1 n2)
			(type bignum-index x y res))
	       (flet ((bigger-p ()
			;; We guarantee in the caller that n1 is either
			;; n2 or (1+ n2)
			(when (> n1 n2)
			  ;;(assert (= n1 (1+ n2)))
			  (if (zerop (%bignum-ref result (+ x n1 -1)))
			      (progn
				(setf (%bignum-ref result (+ res n1 -1)) 0)
				(decf n1))				
			      (return-from bigger-p t)))
			;; Return T if X >= Y.		      
			(do ((i (1- (+ x n1)) (1- i))
			     (j (1- (+ y n2)) (1- j)))
			    (())
			  (declare (type bignum-index i j))
			  (let ((x-digit (%bignum-ref result i))
				(y-digit (%bignum-ref result j)))
			    (when (/= x-digit y-digit)
			      (return-from bigger-p (>= x-digit y-digit))))
			  ;; Need to test for the exit condition here
			  ;; instead of in the DO headers because when
			  ;; X or Y is 0 and the numbers are equal I
			  ;; or J will be decremented to -1, which
			  ;; isn't of type bignum-index.
			  (when (= i x)
			    (return-from bigger-p t)))))
		 (declare (inline bigger-p))
		 (let ((sign 1))
		   (when (not (bigger-p))		     
		     (rotatef x y)
		     (rotatef n1 n2)
		     (setf sign -1))
		   (let ((borrow 1))
		     (do ((i x (1+ i))
			  (j y (1+ j))
			  (k res (1+ k))
			  (end (+ x n1))
			  (end2 (+ y n2)))
			 ((>= i end))
		       (declare (type bignum-index i j k end))
		       (multiple-value-bind (diff new-borrow)
			   (%subtract-with-borrow (%bignum-ref result i)
						  (if (>= j end2)
						      0
						      (%bignum-ref result j))
						  borrow)
			 (setf (%bignum-ref result k) diff)
			 (setf borrow new-borrow)))
		     sign))))
	     (propagate-carry (start end carry)
	       (declare (type bignum-type result)
			(type bignum-index start end)
			(type (integer -4 4) carry))
	       (cond ((minusp carry)
		      (multiple-value-bind (sum borrow)
			  (%subtract-with-borrow (%bignum-ref result start)
						 (- carry) 1)
			(setf (%bignum-ref result start) sum)
			(do ((i (1+ start) (1+ i)))
			    ((or (>= i end) (= borrow 1)) (- borrow 1))
			  (declare (type bignum-index i))
			  (multiple-value-bind (sum next-borrow)
			      (%subtract-with-borrow (%bignum-ref result i)
						     0 borrow)
			    (setf borrow next-borrow)
			    (setf (%bignum-ref result i) sum)))))
		     (t
		      (multiple-value-bind (sum carry)
			  (%add-with-carry (%bignum-ref result start) carry 0)
			(setf (%bignum-ref result start) sum)
			(do ((i (1+ start) (1+ i)))
			    ((or (>= i end) (zerop carry)) carry)
			  (declare (type bignum-index i))
			  (multiple-value-bind (sum next-carry)
			      (%add-with-carry (%bignum-ref result i) 0 carry)
			    (setf carry next-carry)
			    (setf (%bignum-ref result i) sum)))))))
	     (stageC-add ()
	       (let ((carry11 0)
		     (carry12 0)
		     (carry2 0))
		 (do ((i0 ret (1+ i0))
		      (i1 (+ ret d/2) (1+ i1))
		      (i2 (+ ret d) (1+ i2))
		      (i3 (+ ret n/2) (1+ i3))
		      (i4 (+ ret n/2 d/2) (1+ i4))
		      (end (+ ret d/2)))
		     ((>= i0 end))
		   (declare (type bignum-index i0 i1 i2 i3 i4))
		   (let ((b (%bignum-ref result i1)))
		     (multiple-value-bind (sum c1)
			 (%add-with-carry b (%bignum-ref result i0) carry11)
		       (multiple-value-bind (sum c2)
			   (%add-with-carry sum (%bignum-ref result i3) carry12)
			 (setf carry11 c1)
			 (setf carry12 c2)
			 (setf (%bignum-ref result i1) sum)))
		     (multiple-value-bind (sum carry)
			 (%add-with-carry b (%bignum-ref result i4) carry2)
		       (setf (%bignum-ref result i2) sum)
		       (setf carry2 carry))))
		 (+ (propagate-carry (+ ret d) (+ ret n/2) (+ carry11 carry12))
		    carry2)))
	     (stageC-subtract ()
	       (let ((carry11 0)
		     (borrow12 1)
		     (borrow2 1))
		 (do ((i0 ret (1+ i0))
		      (i1 (+ ret d/2) (1+ i1))
		      (i2 (+ ret d) (1+ i2))
		      (i3 (+ ret n/2) (1+ i3))
		      (i4 (+ ret n/2 d/2) (1+ i4))
		      (end (+ ret d/2)))
		     ((>= i0 end))
		   (declare (type bignum-index i0 i1 i2 i3 i4))
		   (let ((b (%bignum-ref result i1)))
		     (multiple-value-bind (sum c1)
			 (%add-with-carry b (%bignum-ref result i0) carry11)
		       (multiple-value-bind (sum b2)
			   (%subtract-with-borrow sum (%bignum-ref result i3)
						  borrow12)
			 (setf carry11 c1)
			 (setf borrow12 b2)
			 (setf (%bignum-ref result i1) sum)))
		     (multiple-value-bind (sum borrow)
			 (%subtract-with-borrow b (%bignum-ref result i4)
						borrow2)
		       (setf (%bignum-ref result i2) sum)
		       (setf borrow2 borrow))))
		 (+ (propagate-carry (+ ret d) (+ ret n/2)
				     (+ carry11 borrow12 -1))
		    (- borrow2 1))))
	     (stageE-add (stageC-carry)
	       ;; Sum the above pieces, result to R0124.
	       (let ((carry1 0)
		     (carry21 0)
		     (carry22 0))
		 (declare (type (mod 2) carry1 carry21 carry22))
		 (do ((i0 ret (1+ i0))
		      (i1 (+ ret d/2) (1+ i1))
		      (i2 (+ ret d) (1+ i2))
		      (i3 (+ ret d d/2) (1+ i3))
		      (i4 (+ ret d d) (1+ i4))
		      (end (+ ret d/2)))
		     ((>= i0 end))
		   (declare (type bignum-index i0 i1 i2 i3 i4))
		   (multiple-value-bind (sum c1)
		       (%add-with-carry (%bignum-ref result i1)
					(%bignum-ref result i3)
					carry1)
		     (setf carry1 c1)
		     (setf (%bignum-ref result i1) sum))
		   (let ((e (%bignum-ref result i4)))
		     (multiple-value-bind (sum c1)
			 (%add-with-carry e (%bignum-ref result i2) carry21)
		       (multiple-value-bind (sum c2)
			   (%add-with-carry sum (%bignum-ref result i3) carry22)
			 (setf carry21 c1)
			 (setf carry22 c2)
			 (setf (%bignum-ref result i2) sum)))
		     (setf (%bignum-ref result i3) e)))
		 (propagate-carry (+ ret d) (+ ret d d) carry1)
		 (propagate-carry (+ ret d d/2)
						 (+ ret d d)
						 (+ carry21 carry22 stageC-carry)))))
      (cond ((<= d (ext:truly-the fixnum *karatsuba-classical-cutoff*))
	     (bignum-multiply-in-place u v ret origd result d))
	    (t
	     (let* ((ret1 (+ ret d/2))
		    (ret3 (+ ret n/2))
		    (u-hi (+ u d/2))
		    (v-hi (+ v d/2))
		    (diff-sign (* (abs-diff d/2 u d/2- u-hi ret result)
				  (abs-diff d/2 v d/2- v-hi ret1 result))))
	       (declare (type bignum-index d/2 ret1 ret3 u-hi v-hi)
			(type (member -1 1) diff-sign))
	       (do ((i ret3 (1+ i)))
		   ((>= i (+ ret3 d)))
		 (setf (%bignum-ref result i) 0))
	       (karatsuba-internal ret ret1 ret3 n/2 d/2 result)
	       (do ((i ret (1+ i)))
		   ((>= i ret3))
		 (setf (%bignum-ref result i) 0))	       
	       (karatsuba-internal u v ret n/2 d/2 result)
	       (let ((stagec-carry (if (minusp diff-sign)
				       (stagec-add)
				       (stagec-subtract))))
		 (setf ret3 (+ ret d d/2))
		 (karatsuba-internal u-hi v-hi ret3 n/2 d/2- result)
		 (when (/= d/2 d/2-)
		   (setf (%bignum-ref result (+ ret3 (* d/2- 2))) 0)
		   (setf (%bignum-ref result (+ ret3 (* d/2- 2) 1)) 0))
		 (stageE-add stageC-carry))))))
    result))


;;; karatsuba  --  Interface.
;;;
;;; Multiply two bignums using the Karatsuba multiplication technique.
;;;
;;; The arguments are currently copied to the end of the working
;;; result causing some extra consing, and the arguments are both
;;; extended to have the same width implementation rather than
;;; exploiting a small width.
;;;
(defun karatsuba (x y)
  (declare (type bignum-type x y)
	   (optimize (speed 3) (safety 0) (debug 3)))
  (flet ((power-of-two (n)
	   ;; Compute the smallest power of two greater than or equal
	   ;; to the given number.
	   (declare (type bignum-index n))
	   (let ((power 1))
	     (declare (type bignum-index power))
	     (loop while (< power n) do
		   (setf power (ash power 1)))
	     power)))
    (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
	   (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
	   (x (if x-plusp x (negate-bignum x)))
	   (y (if y-plusp y (negate-bignum y)))
	   (negate-res (not (eq x-plusp y-plusp)))
	   (x-words (ceiling (integer-length x) 32))
	   (y-words (ceiling (integer-length y) 32))
	   (d (max x-words y-words))
	   (n (power-of-two d)))
      (declare (type bignum-index n))

      (let ((result (%allocate-bignum (* 5 n))))

      ;; Copy X, Y to the end of the result area, so that we can use
      ;; indices off of the result area for accessing everything.
      (do ((k 0 (1+ k))
	   (x-save (* 3 n) (1+ x-save)))
	  ((>= k x-words))
	(setf (%bignum-ref result x-save)
	      (%bignum-ref x k)))
      
      (do ((k 0 (1+ k))
	   (save (* 4 n) (1+ save)))
	  ((>= k y-words))
	(setf (%bignum-ref result save)
	      (%bignum-ref y k)))
      ;; Do it.
      (karatsuba-internal (* 3 n) (* 4 n) 0 (* n 3)
			  (max x-words y-words)
			  result)
      ;; Convert the result into a bignum.
      (let ((new-len (1+ (* 2 n)))
	    (n5 (* n 5)))
	(do ((k (* 2 d) (1+ k)))
	  ;;(assert (< k (%bignum-length result)))
	  (setf (%bignum-ref result k) 0))
	(%bignum-set-length result new-len)
	(when negate-res
	  (negate-bignum-in-place result))
	(%normalize-bignum result (1+ (* 2 n))))))))

(defun classical-multiply-bignums (a b)
wlott's avatar
wlott committed
  (declare (type bignum-type a b))
  (let* ((a-plusp (%bignum-0-or-plusp a (%bignum-length a)))
	 (b-plusp (%bignum-0-or-plusp b (%bignum-length b)))
	 (a (if a-plusp a (negate-bignum a)))
	 (b (if b-plusp b (negate-bignum b)))
	 (len-a (%bignum-length a))
	 (len-b (%bignum-length b))
	 (len-res (+ len-a len-b))
	 (res (%allocate-bignum len-res))
	 (negate-res (not (eq a-plusp b-plusp))))
    (declare (type bignum-index len-a len-b len-res))
wlott's avatar
wlott committed
    (dotimes (i len-a)
      (declare (type bignum-index i))
wlott's avatar
wlott committed
	    (x (%bignum-ref a i))
	    (k i))
	(declare (type bignum-index k)
		 (type bignum-element-type carry-digit x))
	(dotimes (j len-b)
	  (multiple-value-bind (big-carry res-digit)
			       (%multiply-and-add x (%bignum-ref b j)
						  (%bignum-ref res k)
						  carry-digit)
	    (declare (type bignum-element-type big-carry res-digit))
	    (setf (%bignum-ref res k) res-digit)
	    (incf k)))
	(setf (%bignum-ref res k) carry-digit)))
wlott's avatar
wlott committed
    (when negate-res (negate-bignum-in-place res))
    (%normalize-bignum res len-res)))

(defparameter *min-karatsuba-bits* 512
  "Use Karatsuba if the bignums have at least this many bits")

(defun multiply-bignums (a b)
  (declare (type bignum-type a b))
  ;; Use either the classical algorithm or the Karatsuba algorithm
  ;; depending on the size of the numbers.
  (let* ((len-a (integer-length a))
	 (len-b (integer-length b))
	 (min-len (min len-a len-b))
	 (max-len (max len-a len-b)))
    ;; If both numbers are short, use the classical algorithm.  If the
    ;; numbers vary greatly in length, use the classical algorithm.
    ;; Otherwise use the Karatsuba algorithm.
    (if (or (<= min-len *min-karatsuba-bits*)
	    (>= max-len (* 4 min-len)))
	(classical-multiply-bignums a b)
	(karatsuba a b))))

#|
;;;; Test code.

(defun karat-time (n size)
  (declare (fixnum n))
  (let ((max-num (ash 1 size))
	(true 0)
	(kar 0))
    (time
     (dotimes (i n)
       (declare (fixnum i))
       (let ((x (random max-num))
	     (y (random max-num)))
	 (dotimes (j 10)
	   (declare (fixnum j))
	   (setf true (classical-multiply-bignums x y))))))
    (time
     (dotimes (i n)
       (declare (fixnum i))
       (let ((x (random max-num))
	     (y (random max-num)))
	 (dotimes (j 10)
	   (declare (fixnum j))
	   (setf kar (multiply-bignums x y))))))
    (let ((bad (/= kar true)))
      (format t "DIFF~%")
      #+nil
      (when bad
	(format t "(- (* ~X ~X) ~X)~%" x y kar)
	(format t "~X~%~X~%" (* x y) kar)))
    (values)))

(defun karat-time (n size)
  (declare (fixnum n))
  (let* ((max-num (ash 1 size))
	 (true 0)
	 (kar 0)
	 (x (random max-num))
	 (y (random max-num)))
    (ext:gc)
    (format t "~D bits x ~D bits~%" (integer-length x) (integer-length y))
    (time
     (dotimes (i n)
       (declare (fixnum i))
       (dotimes (j 10)
	 (declare (fixnum j))
	 (setf true (classical-multiply-bignums x y)))))
    (ext:gc)
    (time
     (dotimes (i n)
       (declare (fixnum i))
       (dotimes (j 10)
	 (declare (fixnum j))
	 (setf kar (multiply-bignums x y)))))
    (format t "~X~%~X" x y)
    
    (let ((bad (/= kar true)))
      (when bad
	(format t "DIFF~%")
	(format t "(- (* ~X ~X) ~X)~%" x y kar)
	(format t "~X~%~X~%" (* x y) kar)))
    (values)))

(defun karat-time-2 (n size)
  (declare (fixnum n))