Newer
Older
;;; -*- Mode: Lisp; Log: code.log; Package: bignum -*-
;;;
;;; **********************************************************************
;;; 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 $")
;;; **********************************************************************
;;;
;;; This file contains code to implement bignum support.
;;;
(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))
;;;; 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
;;; %BIGNUM-SET-LENGTH
;;; %FIXNUM-DIGIT-WITH-CORRECT-SIGN
;;; %SIGN-DIGIT
;;; %ASHR
;;; %ASHL
;;; %BIGNUM-0-OR-PLUSP
;;; %DIGIT-LOGICAL-SHIFT-RIGHT
;;; 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
;;; Relational operators:
;;; %LOGAND
;;; %LOGIOR
;;; %LOGXOR
;;; LDB
;;; %FIXNUM-TO-DIGIT
;;; TRUNCATE
;;; %FLOOR
;;;
;;;
;;; Note: The floating routines know about the float representation.
;;;
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
;;; 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.
(defconstant maximum-bignum-length (1- (ash 1 (- 32 vm:type-bits))))
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
) ;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)))
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
(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))
;;; %LOGNOT -- Internal.
;;;
(defun %lognot (digit)
(declare (type bignum-element-type digit))
(%lognot digit))
;;; %LOGAND -- Internal.
;;; %LOGIOR -- Internal.
;;; %LOGXOR -- Internal.
;;;
;;; Do the 32bit unsigned op.
;;;
(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))))
#-32x16-divide
;;; %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.
;;;
(let ((a a) (b b) (c c))
(declare (type bignum-element-type a b c))
(setq a a b b c c)
(%floor a b c)))
;;; %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)
(%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))
;;; %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.
;;;
(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.
;;;
(declaim (inline %digit-compare %digit-greater))
(defun %digit-compare (x y)
(= x y))
;;;
(defun %digit-greater (x y)
(> x y))
(declaim (optimize (speed 3) (safety 0)(ext:inhibit-warnings 3)))
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
;;;; 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))
(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))
;;;; 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))
(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))))
(,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)))
(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)
(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))
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
;;; **********************************************************************
;;;
;;; 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)))
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
(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)
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
(%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)
(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))
(let ((carry-digit 0)
(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)
(setf carry-digit big-carry)
(incf k)))
(setf (%bignum-ref res k) carry-digit)))
(when negate-res (negate-bignum-in-place res))
(%normalize-bignum res len-res)))
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
(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))