Loading code/bignum.lisp +135 −131 Original line number Diff line number Diff line Loading @@ -5,7 +5,7 @@ ;;; 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.43 2007/10/10 01:09:32 rtoy Exp $") "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/bignum.lisp,v 1.43.10.1 2008/11/01 16:07:19 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; Loading Loading @@ -830,7 +830,7 @@ down to individual words.") ;;; (defun karatsuba (x y) (declare (type bignum-type x y) (optimize (speed 3) (safety 0) (debug 3))) (optimize (speed 3) (safety 0))) (flet ((power-of-two (n) ;; Compute the smallest power of two greater than or equal ;; to the given number. Loading Loading @@ -3067,74 +3067,7 @@ friends is working. ;;; ;;; These are used by BIGNUM-TRUNCATE and friends in the general case. ;;; (defvar *truncate-x*) (defvar *truncate-y*) ;;; BIGNUM-TRUNCATE -- Public. ;;; ;;; This divides x by y returning the quotient and remainder. In the general ;;; case, we shift y to setup for the algorithm, and we use two buffers to save ;;; consing intermediate values. X gets destructively modified to become the ;;; remainder, and we have to shift it to account for the initial Y shift. ;;; After we multiple bind q and r, we first fix up the signs and then return ;;; the normalized results. ;;; (defun bignum-truncate (x y) (declare (type bignum-type x y)) (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 nil))) (y (if y-plusp y (negate-bignum y nil))) (len-x (%bignum-length x)) (len-y (%bignum-length y))) (multiple-value-bind (q r) (cond ((< len-y 2) (bignum-truncate-single-digit x len-x y)) ((plusp (bignum-compare y x)) (let ((res (%allocate-bignum len-x))) (dotimes (i len-x) (setf (%bignum-ref res i) (%bignum-ref x i))) (values 0 res))) (t (let ((len-x+1 (1+ len-x))) (with-bignum-buffers ((*truncate-x* len-x+1) (*truncate-y* (1+ len-y))) (let ((y-shift (shift-y-for-truncate y))) (shift-and-store-truncate-buffers x len-x y len-y y-shift) (values (do-truncate len-x+1 len-y) ;; DO-TRUNCATE must execute first. (cond ((zerop y-shift) (let ((res (%allocate-bignum len-y))) (declare (type bignum-type res)) (bignum-replace res *truncate-x* :end2 len-y) (%normalize-bignum res len-y))) (t (shift-right-unaligned *truncate-x* 0 y-shift len-y ((= j res-len-1) (setf (%bignum-ref res j) (%ashr (%bignum-ref *truncate-x* i) y-shift)) (%normalize-bignum res res-len)) res))))))))) (let ((quotient (cond ((eq x-plusp y-plusp) q) ((typep q 'fixnum) (the fixnum (- q))) (t (negate-bignum-in-place q)))) (rem (cond (x-plusp r) ((typep r 'fixnum) (the fixnum (- r))) (t (negate-bignum-in-place r))))) (values (if (typep quotient 'fixnum) quotient (%normalize-bignum quotient (%bignum-length quotient))) (if (typep rem 'fixnum) rem (%normalize-bignum rem (%bignum-length rem)))))))) (declaim (ext:start-block bignum-truncate)) ;;; BIGNUM-TRUNCATE-SINGLE-DIGIT -- Internal. ;;; ;;; This divides x by y when y is a single bignum digit. BIGNUM-TRUNCATE fixes Loading Loading @@ -3162,49 +3095,10 @@ friends is working. (setf (%bignum-ref rem 0) r) (values q rem)))) ;;; DO-TRUNCATE -- Internal. ;;; ;;; This divides *truncate-x* by *truncate-y*, and len-x and len-y tell us how ;;; much of the buffers we care about. TRY-BIGNUM-TRUNCATE-GUESS modifies ;;; *truncate-x* on each interation, and this buffer becomes our remainder. ;;; ;;; *truncate-x* definitely has at least three digits, and it has one more than ;;; *truncate-y*. This keeps i, i-1, i-2, and low-x-digit happy. Thanks to ;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS. ;;; (defun do-truncate (len-x len-y) (declare (type bignum-index len-x len-y)) (let* ((len-q (- len-x len-y)) ;; Add one for extra sign digit in case high bit is on. (q (%allocate-bignum (1+ len-q))) (k (1- len-q)) (y1 (%bignum-ref *truncate-y* (1- len-y))) (y2 (%bignum-ref *truncate-y* (- len-y 2))) (i (1- len-x)) (i-1 (1- i)) (i-2 (1- i-1)) (low-x-digit (- i len-y))) (declare (type bignum-index len-q k i i-1 i-2 low-x-digit) (type bignum-element-type y1 y2)) (loop (setf (%bignum-ref q k) (try-bignum-truncate-guess ;; This modifies *truncate-x*. Must access elements each pass. (bignum-truncate-guess y1 y2 (%bignum-ref *truncate-x* i) (%bignum-ref *truncate-x* i-1) (%bignum-ref *truncate-x* i-2)) len-y low-x-digit)) (cond ((zerop k) (return)) (t (decf k) (decf low-x-digit) (shiftf i i-1 i-2 (1- i-2))))) q)) ;;; TRY-BIGNUM-TRUNCATE-GUESS -- Internal. ;;; ;;; This takes a digit guess, multiplies it by *truncate-y* for a result one ;;; greater in length than len-y, and subtracts this result from *truncate-x*. ;;; This takes a digit guess, multiplies it by truncate-y for a result one ;;; greater in length than len-y, and subtracts this result from truncate-x. ;;; Low-x-digit is the first digit of x to start the subtraction, and we know x ;;; is long enough to subtract a len-y plus one length bignum from it. Next we ;;; check the result of the subtraction, and if the high digit in x became Loading @@ -3213,9 +3107,10 @@ friends is working. ;;; subtracting one too many. Knuth shows that the guess is wrong on the order ;;; of 3/b, where b is the base (2 to the digit-size power) -- pretty rarely. ;;; (defun try-bignum-truncate-guess (guess len-y low-x-digit) (defun try-bignum-truncate-guess (guess len-y low-x-digit truncate-x truncate-y) (declare (type bignum-index low-x-digit len-y) (type bignum-element-type guess)) (type bignum-element-type guess) (type bignum-type truncate-x truncate-y)) (let ((carry-digit 0) (borrow 1) (i low-x-digit)) Loading @@ -3225,23 +3120,23 @@ friends is working. ;; Multiply guess and divisor, subtracting from dividend simultaneously. (dotimes (j len-y) (multiple-value-bind (high-digit low-digit) (%multiply-and-add guess (%bignum-ref *truncate-y* j) (%multiply-and-add guess (%bignum-ref truncate-y j) carry-digit) (declare (type bignum-element-type high-digit low-digit)) (setf carry-digit high-digit) (multiple-value-bind (x temp-borrow) (%subtract-with-borrow (%bignum-ref *truncate-x* i) (%subtract-with-borrow (%bignum-ref truncate-x i) low-digit borrow) (declare (type bignum-element-type x) (fixnum temp-borrow)) (setf (%bignum-ref *truncate-x* i) x) (setf (%bignum-ref truncate-x i) x) (setf borrow temp-borrow))) (incf i)) (setf (%bignum-ref *truncate-x* i) (%subtract-with-borrow (%bignum-ref *truncate-x* i) (setf (%bignum-ref truncate-x i) (%subtract-with-borrow (%bignum-ref truncate-x i) carry-digit borrow)) ;; See if guess is off by one, adding one Y back in if necessary. (cond ((%digit-0-or-plusp (%bignum-ref *truncate-x* i)) (cond ((%digit-0-or-plusp (%bignum-ref truncate-x i)) guess) (t ;; If subtraction has negative result, add one divisor value back Loading @@ -3250,17 +3145,58 @@ friends is working. (carry 0)) (dotimes (j len-y) (multiple-value-bind (v k) (%add-with-carry (%bignum-ref *truncate-y* j) (%bignum-ref *truncate-x* i) (%add-with-carry (%bignum-ref truncate-y j) (%bignum-ref truncate-x i) carry) (declare (type bignum-element-type v)) (setf (%bignum-ref *truncate-x* i) v) (setf (%bignum-ref truncate-x i) v) (setf carry k)) (incf i)) (setf (%bignum-ref *truncate-x* i) (%add-with-carry (%bignum-ref *truncate-x* i) 0 carry))) (setf (%bignum-ref truncate-x i) (%add-with-carry (%bignum-ref truncate-x i) 0 carry))) (%subtract-with-borrow guess 1 1))))) ;;; DO-TRUNCATE -- Internal. ;;; ;;; This divides truncate-x by truncate-y, and len-x and len-y tell us how ;;; much of the buffers we care about. TRY-BIGNUM-TRUNCATE-GUESS modifies ;;; truncate-x on each interation, and this buffer becomes our remainder. ;;; ;;; truncate-x definitely has at least three digits, and it has one more than ;;; truncate-y. This keeps i, i-1, i-2, and low-x-digit happy. Thanks to ;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS. ;;; (defun do-truncate (len-x len-y truncate-x truncate-y) (declare (type bignum-index len-x len-y) (type bignum-type truncate-x truncate-y)) (let* ((len-q (- len-x len-y)) ;; Add one for extra sign digit in case high bit is on. (q (%allocate-bignum (1+ len-q))) (k (1- len-q)) (y1 (%bignum-ref truncate-y (1- len-y))) (y2 (%bignum-ref truncate-y (- len-y 2))) (i (1- len-x)) (i-1 (1- i)) (i-2 (1- i-1)) (low-x-digit (- i len-y))) (declare (type bignum-index len-q k i i-1 i-2 low-x-digit) (type bignum-element-type y1 y2)) (loop (setf (%bignum-ref q k) (try-bignum-truncate-guess ;; This modifies truncate-x. Must access elements each pass. (bignum-truncate-guess y1 y2 (%bignum-ref truncate-x i) (%bignum-ref truncate-x i-1) (%bignum-ref truncate-x i-2)) len-y low-x-digit truncate-x truncate-y)) (cond ((zerop k) (return)) (t (decf k) (decf low-x-digit) (shiftf i i-1 i-2 (1- i-2))))) q)) ;;; BIGNUM-TRUNCATE-GUESS -- Internal. ;;; ;;; This returns a guess for the next division step. Y1 is the highest y Loading Loading @@ -3325,6 +3261,7 @@ friends is working. ;;; We shift y to make it sufficiently large that doing the 64-bit by 32-bit ;;; %FLOOR calls ensures the quotient and remainder fit in 32-bits. ;;; (declaim (inline shift-y-for-truncate)) (defun shift-y-for-truncate (y) (let* ((len (%bignum-length y)) (last (%bignum-ref y (1- len)))) Loading @@ -3336,18 +3273,84 @@ friends is working. ;;; ;;; Stores two bignums into the truncation bignum buffers, shifting them on the ;;; way in. This assumes x and y are positive and at least two in length, and ;;; it assumes *truncate-x* and *truncate-y* are one digit longer than x and y. ;;; it assumes truncate-x and truncate-y are one digit longer than x and y. ;;; (defun shift-and-store-truncate-buffers (x len-x y len-y shift) (defun shift-and-store-truncate-buffers (x len-x y len-y shift truncate-x truncate-y) (declare (type bignum-index len-x len-y) (type (integer 0 (#.digit-size)) shift)) (type (integer 0 (#.digit-size)) shift) (type bignum-type truncate-x truncate-y)) (cond ((zerop shift) (bignum-replace *truncate-x* x :end1 len-x) (bignum-replace *truncate-y* y :end1 len-y)) (bignum-replace truncate-x x :end1 len-x) (bignum-replace truncate-y y :end1 len-y)) (t (bignum-ashift-left-unaligned x 0 shift (1+ len-x) truncate-x) (bignum-ashift-left-unaligned y 0 shift (1+ len-y) truncate-y)))) ;;; BIGNUM-TRUNCATE -- Public. ;;; ;;; This divides x by y returning the quotient and remainder. In the general ;;; case, we shift y to setup for the algorithm, and we use two buffers to save ;;; consing intermediate values. X gets destructively modified to become the ;;; remainder, and we have to shift it to account for the initial Y shift. ;;; After we multiple bind q and r, we first fix up the signs and then return ;;; the normalized results. ;;; (defun bignum-truncate (x y) (declare (type bignum-type x y) (optimize (speed 3))) (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 nil))) (y (if y-plusp y (negate-bignum y nil))) (len-x (%bignum-length x)) (len-y (%bignum-length y))) (multiple-value-bind (q r) (cond ((< len-y 2) (bignum-truncate-single-digit x len-x y)) ((plusp (bignum-compare y x)) (let ((res (%allocate-bignum len-x))) (dotimes (i len-x) (setf (%bignum-ref res i) (%bignum-ref x i))) (values 0 res))) (t (let ((len-x+1 (1+ len-x))) (with-bignum-buffers ((truncate-x len-x+1) (truncate-y (1+ len-y))) (let ((y-shift (shift-y-for-truncate y))) (shift-and-store-truncate-buffers x len-x y len-y y-shift truncate-x truncate-y) (values (do-truncate len-x+1 len-y truncate-x truncate-y) ;; DO-TRUNCATE must execute first. (cond ((zerop y-shift) (let ((res (%allocate-bignum len-y))) (declare (type bignum-type res)) (bignum-replace res truncate-x :end2 len-y) (%normalize-bignum res len-y))) (t (bignum-ashift-left-unaligned x 0 shift (1+ len-x) *truncate-x*) (bignum-ashift-left-unaligned y 0 shift (1+ len-y) *truncate-y*)))) (shift-right-unaligned truncate-x 0 y-shift len-y ((= j res-len-1) (setf (%bignum-ref res j) (%ashr (%bignum-ref truncate-x i) y-shift)) (%normalize-bignum res res-len)) res))))))))) (let ((quotient (cond ((eq x-plusp y-plusp) q) ((typep q 'fixnum) (the fixnum (- q))) (t (negate-bignum-in-place q)))) (rem (cond (x-plusp r) ((typep r 'fixnum) (the fixnum (- r))) (t (negate-bignum-in-place r))))) (values (if (typep quotient 'fixnum) quotient (%normalize-bignum quotient (%bignum-length quotient))) (if (typep rem 'fixnum) rem (%normalize-bignum rem (%bignum-length rem)))))))) (declaim (ext:end-block)) ;;;; %FLOOR primitive for BIGNUM-TRUNCATE. Loading @@ -3365,6 +3368,7 @@ friends is working. ;;; %FLOOR for machines with a 32x32 divider. ;;; #+32x16-divide (declaim (inline 32x16-subtract-with-borrow 32x16-add-with-carry 32x16-divide 32x16-multiply 32x16-multiply-split)) Loading code/float.lisp +26 −2 Original line number Diff line number Diff line Loading @@ -5,7 +5,7 @@ ;;; Carnegie Mellon University, and has been placed in the public domain. ;;; (ext:file-comment "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/float.lisp,v 1.41 2008/04/15 16:31:58 rtoy Exp $") "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/float.lisp,v 1.41.8.1 2008/11/01 16:07:20 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; Loading Loading @@ -348,7 +348,31 @@ #+long-float ((long-float) (frob vm:long-float-digits vm:long-float-bias integer-decode-long-denorm))))) integer-decode-long-denorm)) #+double-double ((double-double-float) ;; What exactly is the precision for a double-double? We make ;; it the sum of the precisions of the two components. (let ((hi (double-double-hi f))) (cond ((zerop hi) ;; ANSI CL says the precision is 0 0) ((< (abs hi) (scale-float least-positive-normalized-double-float 53)) ;; For every power of 2 below this, we lose a bit of ;; precision. More or less. (let ((cutoff (nth-value 1 (decode-float (scale-float least-positive-normalized-double-float 53))))) (multiple-value-bind (f exp) (decode-float hi) (declare (ignore f)) (- 106 (- cutoff exp))))) (t ;; Normally we have 106 bits of precision (twice the ;; double-float precision) 106))))))) #+nil Loading code/irrat.lisp +94 −12 Original line number Diff line number Diff line Loading @@ -5,7 +5,7 @@ ;;; Carnegie Mellon University, and has been placed in the public domain. ;;; (ext:file-comment "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.55.8.2 2008/09/26 20:13:10 rtoy Exp $") "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.55.8.2.2.1 2008/11/01 16:07:20 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; Loading Loading @@ -303,6 +303,14 @@ (defparameter *intexp-maximum-exponent* 10000) (define-condition intexp-limit-error (error) ((base :initarg :base :reader intexp-base) (power :initarg :power :reader intexp-power)) (:report (lambda (condition stream) (format stream "The absolute value of ~S exceeds limit ~S." (intexp-power condition) *intexp-maximum-exponent*)))) ;;; This function precisely calculates base raised to an integral power. It ;;; separates the cases by the sign of power, for efficiency reasons, as powers ;;; can be calculated more efficiently if power is a positive integer. Values Loading @@ -316,9 +324,17 @@ (return-from intexp base)) (when (> (abs power) *intexp-maximum-exponent*) (cerror "Continue with calculation." "The absolute value of ~S exceeds ~S." power '*intexp-maximum-exponent* base power)) ;; Allow user the option to continue with calculation, possibly ;; increasing the limit to the given power. (restart-case (error 'intexp-limit-error :base base :power power) (continue () :report "Continue with calculation") (new-limit () :report "Continue with calculation, update limit" (setq *intexp-maximum-exponent* power)))) (cond ((minusp power) (/ (intexp base (- power)))) ((eql base 2) Loading Loading @@ -471,10 +487,19 @@ (coerce (* pow (%cos y*pi)) rtype) (coerce (* pow (%sin y*pi)) rtype))))))))))))) (declare (inline real-expt)) ;; This is really messy and should be cleaned up. The easiest ;; way to see if we're doing what we should is the macroexpand ;; the number-dispatch and check each branch. ;; ;; We try to apply the rule of float precision contagion (CLHS ;; 12.1.4.4): the result has the same precision has the most ;; precise argument. (number-dispatch ((base number) (power number)) (((foreach fixnum (or bignum ratio) (complex rational)) integer) (((foreach fixnum (or bignum ratio) (complex rational)) integer) (intexp base power)) (((foreach single-float double-float) rational) (((foreach single-float double-float) rational) (real-expt base power '(dispatch-type base))) (((foreach fixnum (or bignum ratio) single-float) (foreach ratio single-float)) Loading @@ -485,24 +510,81 @@ ((double-float single-float) (real-expt base power 'double-float)) #+double-double (((foreach fixnum (or bignum ratio) single-float double-float double-double-float) (((foreach fixnum (or bignum ratio) single-float double-float double-double-float) double-double-float) (dd-%pow (coerce base 'double-double-float) power)) #+double-double ((double-double-float (foreach fixnum (or bignum ratio) single-float double-float)) (dd-%pow base (coerce power 'double-double-float))) (((foreach (complex rational) (complex float)) rational) (((foreach (complex rational) (complex single-float) (complex double-float) #+double-double (complex double-double-float)) rational) (* (expt (abs base) power) (cis (* power (phase base))))) (((foreach fixnum (or bignum ratio) single-float double-float #+double-double double-double-float) #+double-double ((double-double-float complex) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (* (log2 base 1w0) (log 2w0)))))) (((foreach fixnum (or bignum ratio) single-float double-float) (foreach (complex double-float))) ;; Result should have double-float accuracy. Use log2 in ;; case the base won't fit in a double-float. (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (* (log2 base) (log 2d0)))))) ((double-float (foreach (complex rational) (complex single-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) #+double-double (((foreach fixnum (or bignum ratio) single-float double-float) (foreach (complex double-double-float))) ;; Result should have double-double-float accuracy. Use log2 ;; in case the base won't fit in a double-float. (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (* (log2 base 1w0) (log 2w0)))))) (((foreach fixnum (or bignum ratio) single-float) (foreach (complex single-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) (((foreach (complex rational) (complex single-float)) (foreach single-float (complex single-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) (((foreach (complex rational) (complex single-float)) (foreach double-float (complex double-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log (coerce base '(complex double-float))))))) #+double-double (((foreach (complex rational) (complex single-float)) (foreach double-double-float (complex double-double-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log (coerce base '(complex double-double-float))))))) (((foreach (complex double-float)) (foreach single-float double-float (complex single-float) (complex double-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) (((foreach (complex float) (complex rational)) (foreach complex double-float single-float #+double-double double-double-float)) #+double-double (((foreach (complex double-float)) (foreach double-double-float (complex double-double-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log (coerce base '(complex double-double-float))))))) #+double-double (((foreach (complex double-double-float)) (foreach float (complex float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))))))) Loading general-info/release-19f.txt +12 −0 Original line number Diff line number Diff line Loading @@ -28,6 +28,8 @@ New in this release: the new network connection o Updated CONNECT-TO-INET-SOCKET to allow binding the newly created socket to a local address - Added UNIX:UNIX-OPENPTY, an interface to the openpty C library function. * ANSI compliance fixes: - Fix bug in backquote printer. If the variable is @foo, we want Loading Loading @@ -87,6 +89,7 @@ New in this release: happen on other architectures. - The interpreter catches invalid EVAL-WHEN situations just like the compiler, instead of silently ignoring them. - FLOAT-PRECISION supports double-double floats. * Trac Tickets: - #16: Read-time hash-table issue Loading @@ -100,6 +103,9 @@ New in this release: fixnums. - #20: Modular arith bug? Workaround applied. - #24: Float contagion for expt Float contagion is applied to the arguments before computing expt. * Other changes: - IS1, IS2, IS3, and IS4 are recognized character names for the Loading @@ -112,6 +118,12 @@ New in this release: - Updated User guide to include more examples of tracing. - Enable gencgc page protection on x86/darwin. This can speed up GC a bit. (Not measured.) - Bignum truncate is significantly faster. Some cl-bench benchmarks are now almost twice as fast. - The continuable error produced by raising an integer to a power exceeding *intexp-maximum-exponent* is now a restart, giving the user the option to continue and update the limit to the new power. * Improvements to the PCL implementation of CLOS: - The compiler and interpreter should handle SLOT-VALUE the same Loading Loading
code/bignum.lisp +135 −131 Original line number Diff line number Diff line Loading @@ -5,7 +5,7 @@ ;;; 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.43 2007/10/10 01:09:32 rtoy Exp $") "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/bignum.lisp,v 1.43.10.1 2008/11/01 16:07:19 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; Loading Loading @@ -830,7 +830,7 @@ down to individual words.") ;;; (defun karatsuba (x y) (declare (type bignum-type x y) (optimize (speed 3) (safety 0) (debug 3))) (optimize (speed 3) (safety 0))) (flet ((power-of-two (n) ;; Compute the smallest power of two greater than or equal ;; to the given number. Loading Loading @@ -3067,74 +3067,7 @@ friends is working. ;;; ;;; These are used by BIGNUM-TRUNCATE and friends in the general case. ;;; (defvar *truncate-x*) (defvar *truncate-y*) ;;; BIGNUM-TRUNCATE -- Public. ;;; ;;; This divides x by y returning the quotient and remainder. In the general ;;; case, we shift y to setup for the algorithm, and we use two buffers to save ;;; consing intermediate values. X gets destructively modified to become the ;;; remainder, and we have to shift it to account for the initial Y shift. ;;; After we multiple bind q and r, we first fix up the signs and then return ;;; the normalized results. ;;; (defun bignum-truncate (x y) (declare (type bignum-type x y)) (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 nil))) (y (if y-plusp y (negate-bignum y nil))) (len-x (%bignum-length x)) (len-y (%bignum-length y))) (multiple-value-bind (q r) (cond ((< len-y 2) (bignum-truncate-single-digit x len-x y)) ((plusp (bignum-compare y x)) (let ((res (%allocate-bignum len-x))) (dotimes (i len-x) (setf (%bignum-ref res i) (%bignum-ref x i))) (values 0 res))) (t (let ((len-x+1 (1+ len-x))) (with-bignum-buffers ((*truncate-x* len-x+1) (*truncate-y* (1+ len-y))) (let ((y-shift (shift-y-for-truncate y))) (shift-and-store-truncate-buffers x len-x y len-y y-shift) (values (do-truncate len-x+1 len-y) ;; DO-TRUNCATE must execute first. (cond ((zerop y-shift) (let ((res (%allocate-bignum len-y))) (declare (type bignum-type res)) (bignum-replace res *truncate-x* :end2 len-y) (%normalize-bignum res len-y))) (t (shift-right-unaligned *truncate-x* 0 y-shift len-y ((= j res-len-1) (setf (%bignum-ref res j) (%ashr (%bignum-ref *truncate-x* i) y-shift)) (%normalize-bignum res res-len)) res))))))))) (let ((quotient (cond ((eq x-plusp y-plusp) q) ((typep q 'fixnum) (the fixnum (- q))) (t (negate-bignum-in-place q)))) (rem (cond (x-plusp r) ((typep r 'fixnum) (the fixnum (- r))) (t (negate-bignum-in-place r))))) (values (if (typep quotient 'fixnum) quotient (%normalize-bignum quotient (%bignum-length quotient))) (if (typep rem 'fixnum) rem (%normalize-bignum rem (%bignum-length rem)))))))) (declaim (ext:start-block bignum-truncate)) ;;; BIGNUM-TRUNCATE-SINGLE-DIGIT -- Internal. ;;; ;;; This divides x by y when y is a single bignum digit. BIGNUM-TRUNCATE fixes Loading Loading @@ -3162,49 +3095,10 @@ friends is working. (setf (%bignum-ref rem 0) r) (values q rem)))) ;;; DO-TRUNCATE -- Internal. ;;; ;;; This divides *truncate-x* by *truncate-y*, and len-x and len-y tell us how ;;; much of the buffers we care about. TRY-BIGNUM-TRUNCATE-GUESS modifies ;;; *truncate-x* on each interation, and this buffer becomes our remainder. ;;; ;;; *truncate-x* definitely has at least three digits, and it has one more than ;;; *truncate-y*. This keeps i, i-1, i-2, and low-x-digit happy. Thanks to ;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS. ;;; (defun do-truncate (len-x len-y) (declare (type bignum-index len-x len-y)) (let* ((len-q (- len-x len-y)) ;; Add one for extra sign digit in case high bit is on. (q (%allocate-bignum (1+ len-q))) (k (1- len-q)) (y1 (%bignum-ref *truncate-y* (1- len-y))) (y2 (%bignum-ref *truncate-y* (- len-y 2))) (i (1- len-x)) (i-1 (1- i)) (i-2 (1- i-1)) (low-x-digit (- i len-y))) (declare (type bignum-index len-q k i i-1 i-2 low-x-digit) (type bignum-element-type y1 y2)) (loop (setf (%bignum-ref q k) (try-bignum-truncate-guess ;; This modifies *truncate-x*. Must access elements each pass. (bignum-truncate-guess y1 y2 (%bignum-ref *truncate-x* i) (%bignum-ref *truncate-x* i-1) (%bignum-ref *truncate-x* i-2)) len-y low-x-digit)) (cond ((zerop k) (return)) (t (decf k) (decf low-x-digit) (shiftf i i-1 i-2 (1- i-2))))) q)) ;;; TRY-BIGNUM-TRUNCATE-GUESS -- Internal. ;;; ;;; This takes a digit guess, multiplies it by *truncate-y* for a result one ;;; greater in length than len-y, and subtracts this result from *truncate-x*. ;;; This takes a digit guess, multiplies it by truncate-y for a result one ;;; greater in length than len-y, and subtracts this result from truncate-x. ;;; Low-x-digit is the first digit of x to start the subtraction, and we know x ;;; is long enough to subtract a len-y plus one length bignum from it. Next we ;;; check the result of the subtraction, and if the high digit in x became Loading @@ -3213,9 +3107,10 @@ friends is working. ;;; subtracting one too many. Knuth shows that the guess is wrong on the order ;;; of 3/b, where b is the base (2 to the digit-size power) -- pretty rarely. ;;; (defun try-bignum-truncate-guess (guess len-y low-x-digit) (defun try-bignum-truncate-guess (guess len-y low-x-digit truncate-x truncate-y) (declare (type bignum-index low-x-digit len-y) (type bignum-element-type guess)) (type bignum-element-type guess) (type bignum-type truncate-x truncate-y)) (let ((carry-digit 0) (borrow 1) (i low-x-digit)) Loading @@ -3225,23 +3120,23 @@ friends is working. ;; Multiply guess and divisor, subtracting from dividend simultaneously. (dotimes (j len-y) (multiple-value-bind (high-digit low-digit) (%multiply-and-add guess (%bignum-ref *truncate-y* j) (%multiply-and-add guess (%bignum-ref truncate-y j) carry-digit) (declare (type bignum-element-type high-digit low-digit)) (setf carry-digit high-digit) (multiple-value-bind (x temp-borrow) (%subtract-with-borrow (%bignum-ref *truncate-x* i) (%subtract-with-borrow (%bignum-ref truncate-x i) low-digit borrow) (declare (type bignum-element-type x) (fixnum temp-borrow)) (setf (%bignum-ref *truncate-x* i) x) (setf (%bignum-ref truncate-x i) x) (setf borrow temp-borrow))) (incf i)) (setf (%bignum-ref *truncate-x* i) (%subtract-with-borrow (%bignum-ref *truncate-x* i) (setf (%bignum-ref truncate-x i) (%subtract-with-borrow (%bignum-ref truncate-x i) carry-digit borrow)) ;; See if guess is off by one, adding one Y back in if necessary. (cond ((%digit-0-or-plusp (%bignum-ref *truncate-x* i)) (cond ((%digit-0-or-plusp (%bignum-ref truncate-x i)) guess) (t ;; If subtraction has negative result, add one divisor value back Loading @@ -3250,17 +3145,58 @@ friends is working. (carry 0)) (dotimes (j len-y) (multiple-value-bind (v k) (%add-with-carry (%bignum-ref *truncate-y* j) (%bignum-ref *truncate-x* i) (%add-with-carry (%bignum-ref truncate-y j) (%bignum-ref truncate-x i) carry) (declare (type bignum-element-type v)) (setf (%bignum-ref *truncate-x* i) v) (setf (%bignum-ref truncate-x i) v) (setf carry k)) (incf i)) (setf (%bignum-ref *truncate-x* i) (%add-with-carry (%bignum-ref *truncate-x* i) 0 carry))) (setf (%bignum-ref truncate-x i) (%add-with-carry (%bignum-ref truncate-x i) 0 carry))) (%subtract-with-borrow guess 1 1))))) ;;; DO-TRUNCATE -- Internal. ;;; ;;; This divides truncate-x by truncate-y, and len-x and len-y tell us how ;;; much of the buffers we care about. TRY-BIGNUM-TRUNCATE-GUESS modifies ;;; truncate-x on each interation, and this buffer becomes our remainder. ;;; ;;; truncate-x definitely has at least three digits, and it has one more than ;;; truncate-y. This keeps i, i-1, i-2, and low-x-digit happy. Thanks to ;;; SHIFT-AND-STORE-TRUNCATE-BUFFERS. ;;; (defun do-truncate (len-x len-y truncate-x truncate-y) (declare (type bignum-index len-x len-y) (type bignum-type truncate-x truncate-y)) (let* ((len-q (- len-x len-y)) ;; Add one for extra sign digit in case high bit is on. (q (%allocate-bignum (1+ len-q))) (k (1- len-q)) (y1 (%bignum-ref truncate-y (1- len-y))) (y2 (%bignum-ref truncate-y (- len-y 2))) (i (1- len-x)) (i-1 (1- i)) (i-2 (1- i-1)) (low-x-digit (- i len-y))) (declare (type bignum-index len-q k i i-1 i-2 low-x-digit) (type bignum-element-type y1 y2)) (loop (setf (%bignum-ref q k) (try-bignum-truncate-guess ;; This modifies truncate-x. Must access elements each pass. (bignum-truncate-guess y1 y2 (%bignum-ref truncate-x i) (%bignum-ref truncate-x i-1) (%bignum-ref truncate-x i-2)) len-y low-x-digit truncate-x truncate-y)) (cond ((zerop k) (return)) (t (decf k) (decf low-x-digit) (shiftf i i-1 i-2 (1- i-2))))) q)) ;;; BIGNUM-TRUNCATE-GUESS -- Internal. ;;; ;;; This returns a guess for the next division step. Y1 is the highest y Loading Loading @@ -3325,6 +3261,7 @@ friends is working. ;;; We shift y to make it sufficiently large that doing the 64-bit by 32-bit ;;; %FLOOR calls ensures the quotient and remainder fit in 32-bits. ;;; (declaim (inline shift-y-for-truncate)) (defun shift-y-for-truncate (y) (let* ((len (%bignum-length y)) (last (%bignum-ref y (1- len)))) Loading @@ -3336,18 +3273,84 @@ friends is working. ;;; ;;; Stores two bignums into the truncation bignum buffers, shifting them on the ;;; way in. This assumes x and y are positive and at least two in length, and ;;; it assumes *truncate-x* and *truncate-y* are one digit longer than x and y. ;;; it assumes truncate-x and truncate-y are one digit longer than x and y. ;;; (defun shift-and-store-truncate-buffers (x len-x y len-y shift) (defun shift-and-store-truncate-buffers (x len-x y len-y shift truncate-x truncate-y) (declare (type bignum-index len-x len-y) (type (integer 0 (#.digit-size)) shift)) (type (integer 0 (#.digit-size)) shift) (type bignum-type truncate-x truncate-y)) (cond ((zerop shift) (bignum-replace *truncate-x* x :end1 len-x) (bignum-replace *truncate-y* y :end1 len-y)) (bignum-replace truncate-x x :end1 len-x) (bignum-replace truncate-y y :end1 len-y)) (t (bignum-ashift-left-unaligned x 0 shift (1+ len-x) truncate-x) (bignum-ashift-left-unaligned y 0 shift (1+ len-y) truncate-y)))) ;;; BIGNUM-TRUNCATE -- Public. ;;; ;;; This divides x by y returning the quotient and remainder. In the general ;;; case, we shift y to setup for the algorithm, and we use two buffers to save ;;; consing intermediate values. X gets destructively modified to become the ;;; remainder, and we have to shift it to account for the initial Y shift. ;;; After we multiple bind q and r, we first fix up the signs and then return ;;; the normalized results. ;;; (defun bignum-truncate (x y) (declare (type bignum-type x y) (optimize (speed 3))) (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 nil))) (y (if y-plusp y (negate-bignum y nil))) (len-x (%bignum-length x)) (len-y (%bignum-length y))) (multiple-value-bind (q r) (cond ((< len-y 2) (bignum-truncate-single-digit x len-x y)) ((plusp (bignum-compare y x)) (let ((res (%allocate-bignum len-x))) (dotimes (i len-x) (setf (%bignum-ref res i) (%bignum-ref x i))) (values 0 res))) (t (let ((len-x+1 (1+ len-x))) (with-bignum-buffers ((truncate-x len-x+1) (truncate-y (1+ len-y))) (let ((y-shift (shift-y-for-truncate y))) (shift-and-store-truncate-buffers x len-x y len-y y-shift truncate-x truncate-y) (values (do-truncate len-x+1 len-y truncate-x truncate-y) ;; DO-TRUNCATE must execute first. (cond ((zerop y-shift) (let ((res (%allocate-bignum len-y))) (declare (type bignum-type res)) (bignum-replace res truncate-x :end2 len-y) (%normalize-bignum res len-y))) (t (bignum-ashift-left-unaligned x 0 shift (1+ len-x) *truncate-x*) (bignum-ashift-left-unaligned y 0 shift (1+ len-y) *truncate-y*)))) (shift-right-unaligned truncate-x 0 y-shift len-y ((= j res-len-1) (setf (%bignum-ref res j) (%ashr (%bignum-ref truncate-x i) y-shift)) (%normalize-bignum res res-len)) res))))))))) (let ((quotient (cond ((eq x-plusp y-plusp) q) ((typep q 'fixnum) (the fixnum (- q))) (t (negate-bignum-in-place q)))) (rem (cond (x-plusp r) ((typep r 'fixnum) (the fixnum (- r))) (t (negate-bignum-in-place r))))) (values (if (typep quotient 'fixnum) quotient (%normalize-bignum quotient (%bignum-length quotient))) (if (typep rem 'fixnum) rem (%normalize-bignum rem (%bignum-length rem)))))))) (declaim (ext:end-block)) ;;;; %FLOOR primitive for BIGNUM-TRUNCATE. Loading @@ -3365,6 +3368,7 @@ friends is working. ;;; %FLOOR for machines with a 32x32 divider. ;;; #+32x16-divide (declaim (inline 32x16-subtract-with-borrow 32x16-add-with-carry 32x16-divide 32x16-multiply 32x16-multiply-split)) Loading
code/float.lisp +26 −2 Original line number Diff line number Diff line Loading @@ -5,7 +5,7 @@ ;;; Carnegie Mellon University, and has been placed in the public domain. ;;; (ext:file-comment "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/float.lisp,v 1.41 2008/04/15 16:31:58 rtoy Exp $") "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/float.lisp,v 1.41.8.1 2008/11/01 16:07:20 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; Loading Loading @@ -348,7 +348,31 @@ #+long-float ((long-float) (frob vm:long-float-digits vm:long-float-bias integer-decode-long-denorm))))) integer-decode-long-denorm)) #+double-double ((double-double-float) ;; What exactly is the precision for a double-double? We make ;; it the sum of the precisions of the two components. (let ((hi (double-double-hi f))) (cond ((zerop hi) ;; ANSI CL says the precision is 0 0) ((< (abs hi) (scale-float least-positive-normalized-double-float 53)) ;; For every power of 2 below this, we lose a bit of ;; precision. More or less. (let ((cutoff (nth-value 1 (decode-float (scale-float least-positive-normalized-double-float 53))))) (multiple-value-bind (f exp) (decode-float hi) (declare (ignore f)) (- 106 (- cutoff exp))))) (t ;; Normally we have 106 bits of precision (twice the ;; double-float precision) 106))))))) #+nil Loading
code/irrat.lisp +94 −12 Original line number Diff line number Diff line Loading @@ -5,7 +5,7 @@ ;;; Carnegie Mellon University, and has been placed in the public domain. ;;; (ext:file-comment "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.55.8.2 2008/09/26 20:13:10 rtoy Exp $") "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.55.8.2.2.1 2008/11/01 16:07:20 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; Loading Loading @@ -303,6 +303,14 @@ (defparameter *intexp-maximum-exponent* 10000) (define-condition intexp-limit-error (error) ((base :initarg :base :reader intexp-base) (power :initarg :power :reader intexp-power)) (:report (lambda (condition stream) (format stream "The absolute value of ~S exceeds limit ~S." (intexp-power condition) *intexp-maximum-exponent*)))) ;;; This function precisely calculates base raised to an integral power. It ;;; separates the cases by the sign of power, for efficiency reasons, as powers ;;; can be calculated more efficiently if power is a positive integer. Values Loading @@ -316,9 +324,17 @@ (return-from intexp base)) (when (> (abs power) *intexp-maximum-exponent*) (cerror "Continue with calculation." "The absolute value of ~S exceeds ~S." power '*intexp-maximum-exponent* base power)) ;; Allow user the option to continue with calculation, possibly ;; increasing the limit to the given power. (restart-case (error 'intexp-limit-error :base base :power power) (continue () :report "Continue with calculation") (new-limit () :report "Continue with calculation, update limit" (setq *intexp-maximum-exponent* power)))) (cond ((minusp power) (/ (intexp base (- power)))) ((eql base 2) Loading Loading @@ -471,10 +487,19 @@ (coerce (* pow (%cos y*pi)) rtype) (coerce (* pow (%sin y*pi)) rtype))))))))))))) (declare (inline real-expt)) ;; This is really messy and should be cleaned up. The easiest ;; way to see if we're doing what we should is the macroexpand ;; the number-dispatch and check each branch. ;; ;; We try to apply the rule of float precision contagion (CLHS ;; 12.1.4.4): the result has the same precision has the most ;; precise argument. (number-dispatch ((base number) (power number)) (((foreach fixnum (or bignum ratio) (complex rational)) integer) (((foreach fixnum (or bignum ratio) (complex rational)) integer) (intexp base power)) (((foreach single-float double-float) rational) (((foreach single-float double-float) rational) (real-expt base power '(dispatch-type base))) (((foreach fixnum (or bignum ratio) single-float) (foreach ratio single-float)) Loading @@ -485,24 +510,81 @@ ((double-float single-float) (real-expt base power 'double-float)) #+double-double (((foreach fixnum (or bignum ratio) single-float double-float double-double-float) (((foreach fixnum (or bignum ratio) single-float double-float double-double-float) double-double-float) (dd-%pow (coerce base 'double-double-float) power)) #+double-double ((double-double-float (foreach fixnum (or bignum ratio) single-float double-float)) (dd-%pow base (coerce power 'double-double-float))) (((foreach (complex rational) (complex float)) rational) (((foreach (complex rational) (complex single-float) (complex double-float) #+double-double (complex double-double-float)) rational) (* (expt (abs base) power) (cis (* power (phase base))))) (((foreach fixnum (or bignum ratio) single-float double-float #+double-double double-double-float) #+double-double ((double-double-float complex) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (* (log2 base 1w0) (log 2w0)))))) (((foreach fixnum (or bignum ratio) single-float double-float) (foreach (complex double-float))) ;; Result should have double-float accuracy. Use log2 in ;; case the base won't fit in a double-float. (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (* (log2 base) (log 2d0)))))) ((double-float (foreach (complex rational) (complex single-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) #+double-double (((foreach fixnum (or bignum ratio) single-float double-float) (foreach (complex double-double-float))) ;; Result should have double-double-float accuracy. Use log2 ;; in case the base won't fit in a double-float. (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (* (log2 base 1w0) (log 2w0)))))) (((foreach fixnum (or bignum ratio) single-float) (foreach (complex single-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) (((foreach (complex rational) (complex single-float)) (foreach single-float (complex single-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) (((foreach (complex rational) (complex single-float)) (foreach double-float (complex double-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log (coerce base '(complex double-float))))))) #+double-double (((foreach (complex rational) (complex single-float)) (foreach double-double-float (complex double-double-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log (coerce base '(complex double-double-float))))))) (((foreach (complex double-float)) (foreach single-float double-float (complex single-float) (complex double-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) (((foreach (complex float) (complex rational)) (foreach complex double-float single-float #+double-double double-double-float)) #+double-double (((foreach (complex double-float)) (foreach double-double-float (complex double-double-float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log (coerce base '(complex double-double-float))))))) #+double-double (((foreach (complex double-double-float)) (foreach float (complex float))) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))))))) Loading
general-info/release-19f.txt +12 −0 Original line number Diff line number Diff line Loading @@ -28,6 +28,8 @@ New in this release: the new network connection o Updated CONNECT-TO-INET-SOCKET to allow binding the newly created socket to a local address - Added UNIX:UNIX-OPENPTY, an interface to the openpty C library function. * ANSI compliance fixes: - Fix bug in backquote printer. If the variable is @foo, we want Loading Loading @@ -87,6 +89,7 @@ New in this release: happen on other architectures. - The interpreter catches invalid EVAL-WHEN situations just like the compiler, instead of silently ignoring them. - FLOAT-PRECISION supports double-double floats. * Trac Tickets: - #16: Read-time hash-table issue Loading @@ -100,6 +103,9 @@ New in this release: fixnums. - #20: Modular arith bug? Workaround applied. - #24: Float contagion for expt Float contagion is applied to the arguments before computing expt. * Other changes: - IS1, IS2, IS3, and IS4 are recognized character names for the Loading @@ -112,6 +118,12 @@ New in this release: - Updated User guide to include more examples of tracing. - Enable gencgc page protection on x86/darwin. This can speed up GC a bit. (Not measured.) - Bignum truncate is significantly faster. Some cl-bench benchmarks are now almost twice as fast. - The continuable error produced by raising an integer to a power exceeding *intexp-maximum-exponent* is now a restart, giving the user the option to continue and update the limit to the new power. * Improvements to the PCL implementation of CLOS: - The compiler and interpreter should handle SLOT-VALUE the same Loading