November 2014 - cmucl-cvs - mailman3.common-lisp.net

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-29-g46b8066
by Raymond Toy 30 Nov '14

30 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 46b8066952521e9d7c43ff101677725425c851a8 (commit) from 70181ee7553a350eedd7acb2815b98dde6ce846c (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 46b8066952521e9d7c43ff101677725425c851a8 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sat Nov 29 21:43:09 2014 -0800 See ticket:91: Forgot to export FLOAT-EXCEPTIONS-SUMMARY-BYTE. diff --git a/src/code/exports.lisp b/src/code/exports.lisp index a6dcc0b..556a20e 100644 --- a/src/code/exports.lisp +++ b/src/code/exports.lisp @@ -1165,6 +1165,7 @@ "COMPATIBLE-FUNCTION-TYPES-P" "FLOAT-INVALID-OP-1-BYTE" "FLOAT-INVALID-OP-2-BYTE" + "FLOAT-EXCEPTIONS-SUMMARY-BYTE" "GET-FP-OPERANDS" "FUSED-MULTIPLY-ADD" "FUSED-MULTIPLY-SUBTRACT") diff --git a/src/compiler/ppc/parms.lisp b/src/compiler/ppc/parms.lisp index cc2c5ba..dda73eb 100644 --- a/src/compiler/ppc/parms.lisp +++ b/src/compiler/ppc/parms.lisp @@ -67,7 +67,8 @@ float-imprecise-trap-bit float-invalid-trap-bit float-divide-by-zero-trap-bit float-invalid-op-1-byte - float-invalid-op-2-byte)) + float-invalid-op-2-byte + float-exceptions-summary-byte)) #+double-double (export '(double-double-float-digits)) ----------------------------------------------------------------------- Summary of changes: src/code/exports.lisp | 1 + src/compiler/ppc/parms.lisp | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-28-g70181ee
by Raymond Toy 30 Nov '14

30 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 70181ee7553a350eedd7acb2815b98dde6ce846c (commit) from 1e4040f32e973a1debcec908e4e69147c31a67ec (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 70181ee7553a350eedd7acb2815b98dde6ce846c Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sat Nov 29 21:00:38 2014 -0800 Partially addresses ticket:91 by allowing the test suite to finish. The tests still fail (along with others), but the test suite will finish. * code/float-trap.lisp: * Update SET-FLOATING-POINT-MODES: * When we clear out bits in float-invalid-op-1-byte, clear out the float-invalid-op-2-byte. Not necessary for this fix, but those bits signal other invalid operations from sqrt and such. * When clearing the exceptions, clear the sticky exceptions summary bit too. * Fix SIGFPE-HANDLER: * As above, clear out the bit in float-invalid-op-2-byte and the exceptions summary bit. * Only clear out the new exception bits, as done for SSE2. * Mask out the float-invalid-op-2-byte and the summary bit before setting the new modes. * compiler/ppc/parms.lisp * Define float-exceptions-summary-byte. * code/exports.lisp: * Export FLOAT-INVALID-OP-2-BYTE diff --git a/src/code/exports.lisp b/src/code/exports.lisp index 99e0762..a6dcc0b 100644 --- a/src/code/exports.lisp +++ b/src/code/exports.lisp @@ -1164,6 +1164,7 @@ (:export "SIGCONTEXT-LR" "COMPATIBLE-FUNCTION-TYPES-P" "FLOAT-INVALID-OP-1-BYTE" + "FLOAT-INVALID-OP-2-BYTE" "GET-FP-OPERANDS" "FUSED-MULTIPLY-ADD" "FUSED-MULTIPLY-SUBTRACT") diff --git a/src/code/float-trap.lisp b/src/code/float-trap.lisp index a736e0d..d67a4ff 100644 --- a/src/code/float-trap.lisp +++ b/src/code/float-trap.lisp @@ -193,7 +193,8 @@ #+(and darwin ppc) (when (member :invalid current-exceptions) ;; Clear out the bits for the detected invalid operation - (setf (ldb vm:float-invalid-op-1-byte modes) 0))) + (setf (ldb vm:float-invalid-op-1-byte modes) 0) + (setf (ldb vm:float-invalid-op-2-byte modes) 0))) (when accrued-x-p (setf (ldb float-sticky-bits modes) @@ -201,12 +202,18 @@ #+(and darwin ppc) (when (member :invalid current-exceptions) ;; Clear out the bits for the detected invalid operation - (setf (ldb vm:float-invalid-op-1-byte modes) 0))) + (setf (ldb vm:float-invalid-op-1-byte modes) 0) + (setf (ldb vm:float-invalid-op-2-byte modes) 0))) + + #+(and darwin ppc) + (when (or accrued-x-p current-x-p) + (setf (ldb vm:float-exceptions-summary-byte modes) 0)) (when fast-mode-p (if fast-mode (setq modes (logior float-fast-bit modes)) (setq modes (logand (lognot float-fast-bit) modes)))) + (setf (floating-point-modes) modes)) (values)) @@ -261,7 +268,12 @@ (traps (logand (ldb float-exceptions-byte modes) (ldb float-traps-byte modes)))) #+(and darwin ppc) - (let ((new-modes modes)) + (let* ((new-modes modes) + (new-exceptions (logandc2 (ldb float-exceptions-byte new-modes) + traps))) + ;; (format t "sigfpe: modes = #B~32,'0b~%" modes) + ;; (format t "sigfpe: new-exc = #B~32,'0b~%" new-exceptions) + (setf (ldb float-exceptions-byte new-modes) new-exceptions) ;; Clear out all exceptions and save them to the context. ;; ;; XXX: Should we just clear out the bits for the traps that are @@ -273,6 +285,9 @@ ;; XXX: Should we only do that if the invalid trap is enabled? (setf (ldb float-invalid-op-1-byte new-modes) 0) (setf (ldb float-invalid-op-2-byte new-modes) 0) + ;; Clear the FP exception summary bit too. + (setf (ldb float-exceptions-summary-byte new-modes) 0) + ;; (format t "sigfpe: new modes = #B~32,'0b~%" new-modes) (setf (floating-point-modes) new-modes) (setf (sigcontext-floating-point-modes (alien:sap-alien scp (* unix:sigcontext))) @@ -371,7 +386,10 @@ ;; an exception. #+ppc (invalid-mask (if (member :invalid traps) - (dpb 0 vm:float-invalid-op-1-byte #xffffffff) + (dpb 0 + (byte 1 31) + (dpb 0 vm::float-invalid-op-2-byte + (dpb 0 vm:float-invalid-op-1-byte #xffffffff))) #xffffffff)) (orig-modes (gensym))) `(let ((,orig-modes (floating-point-modes))) @@ -387,4 +405,4 @@ ,(logand trap-mask exception-mask) #+ppc ,invalid-mask - #+mips ,(dpb 0 float-exceptions-byte #xffffffff)))))))) + #+mips ,(dpb 0 float-exceptions-byte #xffffffff)))))))) diff --git a/src/compiler/ppc/parms.lisp b/src/compiler/ppc/parms.lisp index dab5ce5..cc2c5ba 100644 --- a/src/compiler/ppc/parms.lisp +++ b/src/compiler/ppc/parms.lisp @@ -66,7 +66,8 @@ float-underflow-trap-bit float-overflow-trap-bit float-imprecise-trap-bit float-invalid-trap-bit float-divide-by-zero-trap-bit - float-invalid-op-1-byte)) + float-invalid-op-1-byte + float-invalid-op-2-byte)) #+double-double (export '(double-double-float-digits)) @@ -152,6 +153,10 @@ ;; operation such as some software request, invalid square root, or ;; invalid integer convert. (defconstant float-invalid-op-2-byte (byte 3 8)) +;; FP exception summary bit. This is set if an FP instruction causes +;; any of the FP exception bits to transition from 0 to 1 (meaning +;; that exception occurred.) +(defconstant float-exceptions-summary-byte (byte 1 31)) (defconstant float-fast-bit (ash 1 2)) ; Non-IEEE mode ); eval-when ----------------------------------------------------------------------- Summary of changes: src/code/exports.lisp | 1 + src/code/float-trap.lisp | 28 +++++++++++++++++++++++----- src/compiler/ppc/parms.lisp | 7 ++++++- 3 files changed, 30 insertions(+), 6 deletions(-) hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-27-g1e4040f
by Raymond Toy 29 Nov '14

29 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 1e4040f32e973a1debcec908e4e69147c31a67ec (commit) from 1fc53cc336bf03c52a438bb237d95080f00dd203 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 1e4040f32e973a1debcec908e4e69147c31a67ec Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sat Nov 29 00:12:55 2014 -0800 Add function to print a description of the value of the FPSCR (floating-point status and control register). diff --git a/src/compiler/ppc/float.lisp b/src/compiler/ppc/float.lisp index c5c944b..e928619 100644 --- a/src/compiler/ppc/float.lisp +++ b/src/compiler/ppc/float.lisp @@ -938,6 +938,63 @@ (defknown ((setf floating-point-modes)) (float-modes) float-modes) +(defun describe-floating-point-modes (&optional (modes (floating-point-modes))) + "Print a description of each bit of the floating point modes in MODES." + (format t "31 FX = ~B FP exception summary~%" + (ldb (byte 1 31) modes)) + (format t "30 FEX = ~B FP enabled exception summary~%" + (ldb (byte 1 30) modes)) + (format t "29 VX = ~B FP invalid operation exception summary~%" + (ldb (byte 1 29) modes)) + (format t "28 OX = ~B FP overflow exception (sticky)~%" + (ldb (byte 1 28) modes)) + (format t "27 UX = ~B FP underflow exception (sticky)~%" + (ldb (byte 1 27) modes)) + (format t "26 ZX = ~B FP zero divide exception (sticky)~%" + (ldb (byte 1 26) modes)) + (format t "25 XX = ~B FP inexact exception (sticky)~%" + (ldb (byte 1 25) modes)) + (format t "24 VXSNAN = ~B FP invalid operation with SNaN (sticky)~%" + (ldb (byte 1 24) modes)) + (format t "23 VXISI = ~B FP invalid operation for inf - inf (sticky)~%" + (ldb (byte 1 23) modes)) + (format t "22 VXIDI = ~B FP invalid operation for inf + inf (sticky)~%" + (ldb (byte 1 22) modes)) + (format t "21 VXZDZ = ~B FP invalid operation for 0/0 (sticky)~%" + (ldb (byte 1 21) modes)) + (format t "20 VXIMZ = ~B FP invalid operation for inf * 0 (sticky)~%" + (ldb (byte 1 20) modes)) + (format t "19 VXVC = ~B FP invalid operation for invalid compare (sticky)~%" + (ldb (byte 1 19) modes)) + (format t "18 FR = ~B FP fraction rounded~%" + (ldb (byte 1 18) modes)) + (format t "17 FI = ~B FP fraction inexact~%" + (ldb (byte 1 17) modes)) + (format t "12 FPRF = ~5,'0B FP result flags~%" + (ldb (byte 5 12) modes)) + (format t "11 Reserved~%") + (format t "10 VXSOFT = ~B FP invalid operation for software request (sticky)~%" + (ldb (byte 1 10) modes)) + (format t " 9 VXSQRT = ~B FP invalid operation for sqrt (sticky)~%" + (ldb (byte 1 9) modes)) + (format t " 8 VXCVI = ~B FP invalid operation for integer convert (sticky)~%" + (ldb (byte 1 8) modes)) + (format t " 7 VE = ~B FP invalid operation exception enable~%" + (ldb (byte 1 7) modes)) + (format t " 6 OE = ~B FP overflow exception enable~%" + (ldb (byte 1 6) modes)) + (format t " 5 UE = ~B FP underflow exception enable~%" + (ldb (byte 1 5) modes)) + (format t " 4 ZE = ~B FP zero divide exception enable~%" + (ldb (byte 1 4) modes)) + (format t " 3 XE = ~B FP inexact operation enable~%" + (ldb (byte 1 3) modes)) + (format t " 2 NI = ~B FP non-IEEE mode enable~%" + (ldb (byte 1 2) modes)) + (format t " 0 RN = ~2,'0B FP rounding control + (00 = nearest, 01 = zero, 10 = +inf, 11 = -inf~%" + (ldb (byte 2 0) modes))) + (define-vop (floating-point-modes) (:results (res :scs (unsigned-reg))) (:result-types unsigned-num) ----------------------------------------------------------------------- Summary of changes: src/compiler/ppc/float.lisp | 57 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-26-g1fc53cc
by Raymond Toy 28 Nov '14

28 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 1fc53cc336bf03c52a438bb237d95080f00dd203 (commit) from 2b39eba64aa3f2b122bcc275e6597f9cb2a5f7f4 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 1fc53cc336bf03c52a438bb237d95080f00dd203 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Fri Nov 28 10:40:20 2014 -0800 Signal an error for an invalid :allocation value for a slot. Fixes ticket:95, along with changeset:[2b39eb] diff --git a/src/pcl/std-class.lisp b/src/pcl/std-class.lisp index 1d4049d..3f247f2 100644 --- a/src/pcl/std-class.lisp +++ b/src/pcl/std-class.lisp @@ -435,7 +435,12 @@ slot-name) else do (loop for (option value . more) on slot by #'cddr - if (and (member option '(:allocation :type :initform + if (and (eq option :allocation) + (not (member value '(:class :instance)))) do + (simple-program-error + "~@<Slot ~S: :allocation value must be :class or :instance, not ~S~@:>" + slot-name value) + else if (and (member option '(:type :initform :documentation)) (not (eq unsupplied (getf more option unsupplied)))) do ----------------------------------------------------------------------- Summary of changes: src/pcl/std-class.lisp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-25-g2b39eba
by Raymond Toy 28 Nov '14

28 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 2b39eba64aa3f2b122bcc275e6597f9cb2a5f7f4 (commit) from 55d04a31ba05d20f4e9ab6be9548d374c900b197 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 2b39eba64aa3f2b122bcc275e6597f9cb2a5f7f4 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Fri Nov 28 10:21:52 2014 -0800 Correct the tests for the generic function lambda lists. The check for invalid gf lambda lists happens during macroexpansion, so update the tests to catch this during macroexpansion. This fixes three of the issues in ticket:95. diff --git a/tests/pcl/defgeneric.lisp b/tests/pcl/defgeneric.lisp index 05f01c7..2cf9e07 100644 --- a/tests/pcl/defgeneric.lisp +++ b/tests/pcl/defgeneric.lisp @@ -32,13 +32,12 @@ (in-package "PCL-TESTS") +;; These errors are signaled during macroexpansion of the defgeneric +;; form so catch the errors during macroexpansion. (defmacro define-gf-lambda-list-test (name lambda-list) - `(deftest ,name - (multiple-value-bind (r c) - (ignore-errors - (defgeneric g ,lambda-list)) - (values (null r) (typep c 'error))) - t t)) + `(define-test ,name + (assert-error 'program-error + (macroexpand '(defgeneric g ,lambda-list))))) (define-gf-lambda-list-test defgeneric-lambda-list.0 (a &optional (b 1))) (define-gf-lambda-list-test defgeneric-lambda-list.1 (a &key (b 1))) ----------------------------------------------------------------------- Summary of changes: tests/pcl/defgeneric.lisp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-24-g55d04a3
by Raymond Toy 27 Nov '14

27 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 55d04a31ba05d20f4e9ab6be9548d374c900b197 (commit) from 37d383ae57c2fc9597cb87eac9551af866885cd3 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 55d04a31ba05d20f4e9ab6be9548d374c900b197 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Nov 26 18:35:05 2014 -0800 Fix up a few more issues with log. * src/code/irrat.lisp: * In log10, return the correct type for the result. Previously, we would return a double even if a single were desired. * Fix up a few more cases where we failed to handle log of a huge rational correctly for log base 2 and base 10. * tests/irrat.lisp: * Add some tests to check that log returns the correct type of number for the log base 2 and 10 of very large rationals. diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 325e5bc..0804ef9 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -653,10 +653,10 @@ ;; Number fits in a double, so it's easy (float (%log10 d) float-type)) (t - ;; Number doesn't fit in a double. Do it the hard way. - ;; This should be done more accurately. - (/ (log2 number float-type) - (log2 (float 10 float-type))))))) + ;; Number doesn't fit in a double. Do it the hard way using + ;; log2. This should be done more accurately. + (float (/ (log2 number) + (log2 10)) float-type))))) (defun log (number &optional (base nil base-p)) "Return the logarithm of NUMBER in the base BASE, which defaults to e." @@ -692,10 +692,14 @@ double-float) (log2 number 1d0)) #+double-double - (((foreach integer ratio single-float double-float) + (((foreach single-float double-float) double-double-float) (dd-%log2 (float number 1w0))) #+double-double + (((foreach integer ratio) + double-double-float) + (log2 number 1w0)) + #+double-double ((double-double-float (foreach integer ratio single-float double-float double-double-float)) (dd-%log2 number)))) @@ -722,7 +726,11 @@ double-float) (log10 number 1d0)) #+double-double - (((foreach integer ratio single-float double-float) + (((foreach integer ratio) + ext:double-double-float) + (log10 number 1w0)) + #+double-double + (((foreach single-float double-float) ext:double-double-float) (dd-%log10 (float number 1w0))) #+double-double diff --git a/tests/irrat.lisp b/tests/irrat.lisp index 9538394..a6a405d 100644 --- a/tests/irrat.lisp +++ b/tests/irrat.lisp @@ -88,7 +88,12 @@ (let ((y (kernel::dd-%log2 (sqrt 0.5w0)))) (assert-true (<= (relerr y -1/2) (* 2.7 (scale-float 1d0 (- (float-digits 1w0))))) - y))) + y)) + (assert-true (typep (log (ash 1 3000) 2) 'single-float)) + (assert-true (typep (log (ash 1 3000) 2f0) 'single-float)) + (assert-true (typep (log (ash 1 3000) 2d0) 'double-float)) + (assert-true (typep (log (ash 1 3000) 2w0) 'ext:double-double-float))) + (define-test dd-log2.powers-of-2 (loop for k from -1074 below 1024 @@ -100,7 +105,11 @@ (define-test dd-log10.special-cases (let ((y (kernel::dd-%log10 (sqrt 10w0)))) (assert-true (<= (relerr y 1/2) - (* 0.25 (scale-float 1d0 (- (float-digits 1w0)))))))) + (* 0.25 (scale-float 1d0 (- (float-digits 1w0))))))) + (assert-true (typep (log (ash 1 3000) 10) 'single-float)) + (assert-true (typep (log (ash 1 3000) 10f0) 'single-float)) + (assert-true (typep (log (ash 1 3000) 10d0) 'double-float)) + (assert-true (typep (log (ash 1 3000) 10w0) 'ext:double-double-float))) (define-test dd-log10.powers-of-ten ;; It would be nice if dd-%log10 produce the exact result for powers ----------------------------------------------------------------------- Summary of changes: src/code/irrat.lisp | 20 ++++++++++++++------ tests/irrat.lisp | 13 +++++++++++-- 2 files changed, 25 insertions(+), 8 deletions(-) hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-23-g37d383a
by Raymond Toy 26 Nov '14

26 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 37d383ae57c2fc9597cb87eac9551af866885cd3 (commit) from 0a522bd0c8c38f6ee76f0cc2122d9984c4e269aa (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 37d383ae57c2fc9597cb87eac9551af866885cd3 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Nov 26 14:04:29 2014 -0800 Fix some issues with log2 and log10 of big rationals. * src/code/irrat.lisp: * Add log10 function to handle the computation of the base 10 log of big rationals that might not fit in a double-float. * Fix some issues where CL:LOG wasn't handling logs of big rationals. (A regression). * tests/trac.lisp: * Update trac.8 test to include logs base 10. diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 568fa46..325e5bc 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -645,6 +645,19 @@ (log2-bignum x) (+ n frac)))))) +;; Handle the case where number could be so large that it doesn't fit +;; in a double-float. +(defun log10 (number &optional (float-type 1d0)) + (let ((d (ignore-errors (float number 1d0)))) + (cond (d + ;; Number fits in a double, so it's easy + (float (%log10 d) float-type)) + (t + ;; Number doesn't fit in a double. Do it the hard way. + ;; This should be done more accurately. + (/ (log2 number float-type) + (log2 (float 10 float-type))))))) + (defun log (number &optional (base nil base-p)) "Return the logarithm of NUMBER in the base BASE, which defaults to e." (if base-p @@ -667,12 +680,17 @@ ((double-float (foreach integer ratio single-float double-float)) (%log2 number)) - (((foreach integer ratio single-float) + ((single-float double-float) + (%log2 (float number 1d0))) + ((single-float (foreach integer ratio single-float)) (float (%log2 (float number 1d0)) 1f0)) - (((foreach integer ratio single-float) + (((foreach integer ratio) + (foreach integer ratio single-float)) + (log2 number 1f0)) + (((foreach integer ratio) double-float) - (%log2 (float number 1d0))) + (log2 number 1d0)) #+double-double (((foreach integer ratio single-float double-float) double-double-float) @@ -690,14 +708,19 @@ ;; log base 2 and 10 for non-negative arguments. (number-dispatch ((number real) (base real)) ((double-float - (foreach rational single-float double-float)) + (foreach integer ratio single-float double-float)) (%log10 number)) - (((foreach integer ratio single-float) + ((single-float double-float) + (%log10 (float number 1d0))) + ((single-float (foreach integer ratio single-float)) (float (%log10 (float number 1d0)) 1f0)) - (((foreach integer ratio single-float) + (((foreach integer ratio) + (foreach integer ratio single-float)) + (log10 number 1f0)) + (((foreach integer ratio) double-float) - (%log10 (float number 1d0))) + (log10 number 1d0)) #+double-double (((foreach integer ratio single-float double-float) ext:double-double-float) diff --git a/tests/trac.lisp b/tests/trac.lisp index 2bf420d..198a24d 100644 --- a/tests/trac.lisp +++ b/tests/trac.lisp @@ -22,7 +22,7 @@ (:tag :trac) (assert-false (let (failures) - (dolist (base (list nil 2 2.0 2d0 + (dolist (base (list nil 2 2.0 2d0 10 10.0 10d0 (ash 1 99) (ash 1 3000) 8/7 (/ (ash 1 3000) 7))) (dolist (number (list 100 100.0 100d0 ----------------------------------------------------------------------- Summary of changes: src/code/irrat.lisp | 37 ++++++++++++++++++++++++++++++------- tests/trac.lisp | 2 +- 2 files changed, 31 insertions(+), 8 deletions(-) hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-22-g0a522bd
by Raymond Toy 26 Nov '14

26 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 0a522bd0c8c38f6ee76f0cc2122d9984c4e269aa (commit) via cc8c049fe257b28f470baafcf92b6ceb929582de (commit) from 9918ab2d5794ac01efe17b808b351e25519dc88a (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 0a522bd0c8c38f6ee76f0cc2122d9984c4e269aa Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Nov 26 10:24:17 2014 -0800 Use new log2 function in C * code/exports.lisp: * Export %LOG2. * code/irrat.lisp: * Define %log2 * Use %log2 instead of log2. (This needs work) * compiler/float-tran.lisp: * Use %log2 instead of log2 in the deftransforms. diff --git a/src/code/exports.lisp b/src/code/exports.lisp index e75e5d7..99e0762 100644 --- a/src/code/exports.lisp +++ b/src/code/exports.lisp @@ -2095,7 +2095,7 @@ "%CALLER-FRAME-AND-PC" "%CHECK-BOUND" "%CLOSURE-FUNCTION" "%CLOSURE-INDEX-REF" "%COS" "%COSH" "%DEPOSIT-FIELD" "%DOUBLE-FLOAT" "%DPB" "%EXP" "%EXPM1" "%HYPOT" "%LDB" - "%LOG" "%LOGB" "%LOG10" "%LOG1P" "%LONG-FLOAT" + "%LOG" "%LOGB" "%LOG10" "%LOG1P" "%LOG2" "%LONG-FLOAT" "%MAKE-COMPLEX" "%MAKE-FUNCALLABLE-INSTANCE" "%MAKE-RATIO" "%MASK-FIELD" "%NEGATE" "%POW" "%RAW-BITS" "%RAW-REF-COMPLEX-DOUBLE" "%RAW-REF-COMPLEX-LONG" diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 1179515..568fa46 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -82,6 +82,7 @@ (def-math-rtn ("__ieee754_exp" %exp) 1) (def-math-rtn ("__ieee754_log" %log) 1) (def-math-rtn ("__ieee754_log10" %log10) 1) +(def-math-rtn ("cmucl_log2" %log2) 1) (def-math-rtn ("__ieee754_pow" %pow) 2) #-(or x86 sparc-v7 sparc-v8 sparc-v9) @@ -665,13 +666,13 @@ (number-dispatch ((number real) (base real)) ((double-float (foreach integer ratio single-float double-float)) - (log2 number)) + (%log2 number)) (((foreach integer ratio single-float) (foreach integer ratio single-float)) - (float (log2 (float number 1d0)) 1f0)) + (float (%log2 (float number 1d0)) 1f0)) (((foreach integer ratio single-float) double-float) - (log2 (float number 1d0))) + (%log2 (float number 1d0))) #+double-double (((foreach integer ratio single-float double-float) double-double-float) diff --git a/src/compiler/float-tran.lisp b/src/compiler/float-tran.lisp index 6efe701..c30e3f7 100644 --- a/src/compiler/float-tran.lisp +++ b/src/compiler/float-tran.lisp @@ -865,7 +865,7 @@ (cond ((= y-val 10) `(coerce (kernel:%log10 (float x 1d0)) 'single-float)) ((= y-val 2) - `(coerce (kernel::log2 (float x 1d0)) 'single-float))))) + `(coerce (kernel:%log2 (float x 1d0)) 'single-float))))) (deftransform log ((x y) ((or (member 0d0) (double-float 0d0)) (constant-argument number)) @@ -879,7 +879,7 @@ (cond ((= y-val 10) `(kernel:%log10 (float x 1d0))) ((= y-val 2) - `(kernel::log2 (float x 1d0)))))) + `(kernel:%log2 (float x 1d0)))))) ;;; Handle some simple transformations commit cc8c049fe257b28f470baafcf92b6ceb929582de Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Nov 26 10:21:50 2014 -0800 Extract the log2 code out of e_pow.c to implement log2 function. * log2.c: * New file containing the parts of e_pow.c that implemented a log2 function for use in pow(). * GNUmakefile: * Compile log2.c as part of the build. diff --git a/src/lisp/GNUmakefile b/src/lisp/GNUmakefile index d8fb690..f67a20e 100644 --- a/src/lisp/GNUmakefile +++ b/src/lisp/GNUmakefile @@ -13,7 +13,8 @@ FDLIBM = k_sin.c k_cos.c k_tan.c s_sin.c s_cos.c s_tan.c sincos.c \ e_atan2.c \ e_rem_pio2.c k_rem_pio2.c \ e_log10.c s_scalbn.c \ - setexception.c + setexception.c \ + log2.c SRCS = lisp.c coreparse.c alloc.c monitor.c print.c interr.c \ vars.c parse.c interrupt.c search.c validate.c globals.c \ diff --git a/src/lisp/log2.c b/src/lisp/log2.c new file mode 100644 index 0000000..d56bc43 --- /dev/null +++ b/src/lisp/log2.c @@ -0,0 +1,139 @@ +/* + * ==================================================== + * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. + * + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * This is the log2 algorithm pulled out of __ieee754_pow in e_pow.c + */ + +/* + * Method: Let x = 2 * (1+f) + * 1. Compute and return log2(x) in two pieces: + * log2(x) = w1 + w2, + * where w1 has 53-24 = 29 bit trailing zeros. + */ + +#include "fdlibm.h" + +static double +bp[] = {1.0, 1.5,}, +dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ +dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ +zero = 0.0, +one = 1.0, +two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ + /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ +L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ +L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ +L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ +L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ +L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ +L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ +cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ +cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ +cp_l = -7.02846165095275826516e-09; /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ + +double cmucl_log2(double x) +{ + double ax; + int k, hx, lx, ix; + union { int i[2]; double d; } ux; + + ux.d = x; + hx = ux.i[HIWORD]; lx = ux.i[LOWORD]; + ix = hx&0x7fffffff; + + ax = fabs(x); + + /* special value of x */ + + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx&0x7fffffff)|lx)==0) { + /* log(+-0)=-inf */ + return fdlibm_setexception(-1.0, FDLIBM_DIVIDE_BY_ZERO); + } + + if (hx<0) { + /* log(-#) = NaN */ + return fdlibm_setexception(x, FDLIBM_INVALID); + } + } + + { + double ss,s2,s_h,s_l,t_h,t_l; + double z_h,z_l,p_h,p_l; + double r, u, v, t, t1, t2; + int n, j; + + n = 0; + /* take care subnormal number */ + if(ix<0x00100000) { + ax *= two53; + n -= 53; + ux.d = ax; + ix = ux.i[HIWORD]; + } + n += ((ix)>>20)-0x3ff; + j = ix&0x000fffff; + /* determine interval */ + ix = j|0x3ff00000; /* normalize ix */ + if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ + else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ + else {k=0;n+=1;ix -= 0x00100000;} + ux.d = ax; + ux.i[HIWORD] = ix; + ax = ux.d; + + /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ + u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ + v = one/(ax+bp[k]); + ss = u*v; + s_h = ss; + ux.d = s_h; + ux.i[LOWORD] = 0; + s_h = ux.d; + /* t_h=ax+bp[k] High */ + t_h = zero; + ux.d = t_h; + ux.i[HIWORD]=((ix>>1)|0x20000000)+0x00080000+(k<<18); + t_h = ux.d; + t_l = ax - (t_h-bp[k]); + s_l = v*((u-s_h*t_h)-s_h*t_l); + /* compute log(ax) */ + s2 = ss*ss; + r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); + r += s_l*(s_h+ss); + s2 = s_h*s_h; + t_h = 3.0+s2+r; + ux.d = t_h; + ux.i[LOWORD] = 0; + t_h = ux.d; + t_l = r-((t_h-3.0)-s2); + /* u+v = ss*(1+...) */ + u = s_h*t_h; + v = s_l*t_h+t_l*ss; + /* 2/(3log2)*(ss+...) */ + p_h = u+v; + ux.d = p_h; + ux.i[LOWORD] = 0; + p_h = ux.d; + p_l = v-(p_h-u); + z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ + z_l = cp_l*p_h+p_l*cp+dp_l[k]; + /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ + t = (double)n; + t1 = (((z_h+z_l)+dp_h[k])+t); + ux.d = t1; + ux.i[LOWORD] = 0; + t1 = ux.d; + t2 = z_l-(((t1-t)-dp_h[k])-z_h); + + return t1 + t2; + } +} ----------------------------------------------------------------------- Summary of changes: src/code/exports.lisp | 2 +- src/code/irrat.lisp | 7 ++- src/compiler/float-tran.lisp | 4 +- src/lisp/GNUmakefile | 3 +- src/lisp/log2.c | 139 +++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 148 insertions(+), 7 deletions(-) create mode 100644 src/lisp/log2.c hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-20-g9918ab2
by Raymond Toy 25 Nov '14

25 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 9918ab2d5794ac01efe17b808b351e25519dc88a (commit) via d46a4bfeff75f685f2ccc4a2627a921e46547c1c (commit) via 8e0c67d0c74e1dd5206d2b068734d863440ca286 (commit) from 4d3255aa1a770f59d2851fd2c85707164ca485f5 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 9918ab2d5794ac01efe17b808b351e25519dc88a Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Nov 25 09:35:38 2014 -0800 * Add tests for dd-%log2 and dd-%log10. * Fix the log10.result-types test because we return correctly rounded results for these few tests. diff --git a/tests/irrat.lisp b/tests/irrat.lisp index f1913d3..9538394 100644 --- a/tests/irrat.lisp +++ b/tests/irrat.lisp @@ -5,9 +5,13 @@ (in-package "IRRAT-TESTS") +(defun relerr (actual expected) + (/ (abs (- actual expected)) + expected)) + ;; This tests that log base 2 returns the correct value and the ;; correct type. -(define-test log2 +(define-test log2.result-types (dolist (number '(4 4f0 4d0 #+double-double 4w0)) (dolist (base '(2 2f0 2d0 #+double-double 2w0)) ;; This tests that log returns the correct value and the correct type. @@ -35,54 +39,84 @@ ;; This tests that log base 10 returns the correct value and the ;; correct type. -(define-test log10 +(define-test log10.result-types (dolist (number '(100 100f0 100d0 #+double-double 100w0)) (dolist (base '(10 10f0 10d0 #+double-double 10w0)) ;; This tests that log returns the correct value and the correct type. (let* ((result (log number base)) - (relerr (/ (abs (- result 2)) 2))) - ;; Figure out the expected type of the result and the maximum - ;; allowed relative error. It turns out that for these test - ;; cases, the result is exactly 2 except when the result type - ;; is a double-double-float. In that case, there is a slight - ;; error for (log 100w0 10). - (multiple-value-bind (true-type allowed-error) - (etypecase number - ((or integer single-float) - (etypecase base + (true-type + (etypecase number ((or integer single-float) - (values 'single-float 0)) + (etypecase base + ((or integer single-float) + 'single-float) + (double-float + 'double-float) + #+double-double + (ext:double-double-float + 'ext:double-double-float))) (double-float - (values 'double-float 0)) - #+double-double - (ext:double-double-float - (values 'ext:double-double-float - 7.5d-33)))) - (double-float - (etypecase base - ((or integer single-float double-float) - (values 'double-float 0)) + (etypecase base + ((or integer single-float double-float) + 'double-float) + #+double-double + (ext:double-double-float + 'ext:double-double-float))) #+double-double (ext:double-double-float - (values 'ext:double-double-float - 7.5d-33)))) - #+double-double - (ext:double-double-float - (values 'ext:double-double-float - 7.5d-33))) - (assert-true (<= relerr allowed-error) - number base result relerr allowed-error) - (assert-true (typep result true-type) - number baes result true-type)))))) + 'ext:double-double-float)))) + (assert-equalp 2 result + number base result) + (assert-true (typep result true-type) + number base result true-type))))) -(define-test dd-log2 +(define-test dd-log2.special-cases ;; Verify that for x = 10^k for k = 1 to 300 that (kernel::dd-%log2 ;; x) is close to the expected value. Previously, a bug caused ;; (kernel::dd-%log2 100w0) to give 6.1699... instead of 6.64385. (loop for k from 1 below 300 - and x = (expt 10 k) - and y = (kernel::dd-%log2 (float x 1w0)) - and z = (/ (log (float x 1d0)) (log 2d0)) - and e = (/ (abs (- y z)) z) + for x = (expt 10 k) + for y = (kernel::dd-%log2 (float x 1w0)) + for z = (/ (log (float x 1d0)) (log 2d0)) + for e = (/ (abs (- y z)) z) do (assert-true (<= e 2d-16) - k y z e))) \ No newline at end of file + k y z e)) + (let ((y (kernel::dd-%log2 (sqrt 2w0)))) + (assert-true (<= (relerr y 1/2) + (* 2.7 (scale-float 1d0 (- (float-digits 1w0))))) + y)) + (let ((y (kernel::dd-%log2 (sqrt 0.5w0)))) + (assert-true (<= (relerr y -1/2) + (* 2.7 (scale-float 1d0 (- (float-digits 1w0))))) + y))) + +(define-test dd-log2.powers-of-2 + (loop for k from -1074 below 1024 + for x = (scale-float 1w0 k) + for y = (kernel::dd-%log2 x) + do (assert-equalp k y + k x y))) + +(define-test dd-log10.special-cases + (let ((y (kernel::dd-%log10 (sqrt 10w0)))) + (assert-true (<= (relerr y 1/2) + (* 0.25 (scale-float 1d0 (- (float-digits 1w0)))))))) + +(define-test dd-log10.powers-of-ten + ;; It would be nice if dd-%log10 produce the exact result for powers + ;; of ten, but we currently don't. But note that the maximum + ;; relative error is less than a double-double epsilon. + (let ((threshold (* 0.109 (scale-float 1d0 (- (float-digits 1w0)))))) + (loop for k from -323 below 0 + for x = (expt 10 k) + for y = (kernel::dd-%log10 (float x 1w0)) + for e = (relerr y k) + do (assert-true (<= e threshold) + k e x y)) + (loop for k from 1 to 308 + for x = (expt 10 k) + for y = (kernel::dd-%log10 (float x 1w0)) + for e = (relerr y k) + do (assert-true (<= e threshold) + k e x y)))) + commit d46a4bfeff75f685f2ccc4a2627a921e46547c1c Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Nov 25 09:33:58 2014 -0800 Use log2 and log10 functions when possible instead of using the general case. diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index fa89142..1179515 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -653,31 +653,35 @@ ((and (realp number) (realp base)) (cond ((and (= base 2) - (floatp number) - #+double-double - (not (typep number 'ext:double-double-float)) + ;;(floatp number) (or (plusp number) (eql number 0.0) - (eql number 0d0))) - ;; Do the same thing as the deftranform does for - ;; log base 2 and 10 for non-negative arguments. + (eql number 0d0) + #+double-double + (eql number 0w0))) + ;; Do the same thing as the deftranform does for log + ;; base 2 and 10 for non-negative arguments: handle + ;; the case where number > 0 or equal to +0. (number-dispatch ((number real) (base real)) ((double-float - (foreach integer single-float double-float)) + (foreach integer ratio single-float double-float)) (log2 number)) - ((single-float - (foreach integer single-float)) + (((foreach integer ratio single-float) + (foreach integer ratio single-float)) (float (log2 (float number 1d0)) 1f0)) - ((single-float double-float) + (((foreach integer ratio single-float) + double-float) (log2 (float number 1d0))) #+double-double - (((foreach integer single-float double-float) - ext:double-double-float) - (log2 (float number 1w0) base)))) + (((foreach integer ratio single-float double-float) + double-double-float) + (dd-%log2 (float number 1w0))) + #+double-double + ((double-double-float + (foreach integer ratio single-float double-float double-double-float)) + (dd-%log2 number)))) ((and (= base 10) - (floatp number) - #+double-double - (not (typep number 'double-double-float)) + ;;(floatp number) (or (plusp number) (eql number 0.0) (eql number 0d0))) @@ -685,19 +689,22 @@ ;; log base 2 and 10 for non-negative arguments. (number-dispatch ((number real) (base real)) ((double-float - (foreach double-float single-float integer)) + (foreach rational single-float double-float)) (%log10 number)) - ((single-float - (foreach single-float integer)) + (((foreach integer ratio single-float) + (foreach integer ratio single-float)) (float (%log10 (float number 1d0)) 1f0)) - ((single-float double-float) + (((foreach integer ratio single-float) + double-float) (%log10 (float number 1d0))) #+double-double - (((foreach integer single-float double-float) + (((foreach integer ratio single-float double-float) ext:double-double-float) - ;; This could be more accurate! - (/ (log (float number 1w0)) - (log 10w0))))) + (dd-%log10 (float number 1w0))) + #+double-double + ((double-double-float + (foreach integer ratio single-float double-float double-double-float)) + (dd-%log10 number)))) (t ;; CLHS 12.1.4.1 says ;; commit 8e0c67d0c74e1dd5206d2b068734d863440ca286 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Nov 25 09:32:32 2014 -0800 Return -infinity for %log2 and %log10 of +0. diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp index d0acdd7..36051fc 100644 --- a/src/code/irrat-dd.lisp +++ b/src/code/irrat-dd.lisp @@ -1379,6 +1379,9 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 (declare (type double-double-float x) (optimize (speed 3) (space 0) (inhibit-warnings 3))) + (when (eql x 0w0) + ;; log2(+0) = -infinity + (return-from dd-%log2 (/ -1 x))) (multiple-value-bind (e x y z) (compute-log x) ;; Multiply log of fraction by log2(e) and base 2 exponent by 1 @@ -1395,6 +1398,9 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 (declare (type double-double-float x) (optimize (speed 3) (space 0) (inhibit-warnings 3))) + (when (eql x 0w0) + ;; log2(+0) = -infinity + (return-from dd-%log10 (/ -1 x))) (multiple-value-bind (e x y z) (compute-log x) ;; Multiply log of fraction by log10(e) and base 2 exponent by log10(2). ----------------------------------------------------------------------- Summary of changes: src/code/irrat-dd.lisp | 6 +++ src/code/irrat.lisp | 55 ++++++++++++++----------- tests/irrat.lisp | 110 ++++++++++++++++++++++++++++++++----------------- 3 files changed, 109 insertions(+), 62 deletions(-) hooks/post-receive -- CMU Common Lisp

1 0

[Cmucl-cvs] [git] CMU Common Lisp branch master updated. snapshot-2014-11-17-g4d3255a
by Raymond Toy 24 Nov '14

24 Nov '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, master has been updated via 4d3255aa1a770f59d2851fd2c85707164ca485f5 (commit) from 3f063954c98d21ea8a95388d01db96a1e056c34d (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 4d3255aa1a770f59d2851fd2c85707164ca485f5 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Nov 24 15:20:12 2014 -0800 Add log10 implementation for double-doubles. Since log2 and log10 use basically the same natural log implementation, factor that out the common part into its own routine. diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp index 45b422e..d0acdd7 100644 --- a/src/code/irrat-dd.lisp +++ b/src/code/irrat-dd.lisp @@ -44,6 +44,20 @@ 4.4269504088896340735992468100189213742664595w-1 _N"log2(e)-1") +;; l102a+log102b = log10(2) to extra precision +(defconstant l102a + 0.3125w0) + +(defconstant l102b + -1.14700043360188047862611052755069732318101185w-2) + +;; l10ea + l10eb = log10(2) to extra precsion +(defconstant l10ea + 0.5w0) + +(defconstant l10eb + -6.570551809674817234887108108339491770560299w-2) + (defconstant dd-pi 3.141592653589793238462643383279502884197169w0 _N"Pi") @@ -1241,8 +1255,8 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 (values (- (dd-%%cos reduced)) (dd-%%sin reduced)))))))) -;;; dd-%log2 -;;; Base 2 logarithm. +;;; dd-%log2 and dd-%log10 +;;; Base 2 and base 10 logarithm. ;;; ;;; The argument is separated into its exponent and fractional ;;; parts. If the exponent is between -1 and +1, the (natural) @@ -1254,6 +1268,9 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 ;;; ;;; log(x) = z + z**3 R(z)/S(z). ;;; +;;; This gives the natural log. To get the base 2 and base 10 log, +;;; carefully multiply the natural log by log2(e) or log10(e) as +;;; appropriate. (let ((P (make-array 13 :element-type 'double-double-float :initial-contents '( @@ -1314,52 +1331,82 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 -1.332535117259762928288745111081235577029w6 1.701761051846631278975701529965589676574w6 )))) + (flet + ((compute-log (x) + ;; Common routine to extract the exponent and fraction from x + ;; and compute the log(f). Return the exponent, the fraction, + ;; and the parts of the polynomial computation that is needed + ;; to finish off the final logarithm. + (declare (type double-double-float x) + (optimize (speed 3) (space 0) + (inhibit-warnings 3))) + (multiple-value-bind (x e) + (decode-float x) + (declare (type double-double-float x) + (type double-float-exponent e)) + (let ((z 0w0) + (y 0w0)) + (declare (type double-double-float z y)) + (cond ((or (> e 2) + (< e -2)) + (cond ((< x sqrt-1/2) + ;; 2*(2*x-1)/(2*x+1) + (decf e) + (setf z (- x 0.5w0)) + (setf y (+ (* 0.5w0 z) 0.5w0))) + (t + ;; 2*(x-1)/(x+1) + (setf z (- x 0.5w0)) + (decf z 0.5w0) + (setf y (+ (* 0.5w0 x) 0.5w0)))) + (setf x (/ z y)) + (setf z (* x x)) + (setf y (* x (/ (* z (poly-eval z r)) + (poly-eval-1 z s))))) + (t + (cond ((< x sqrt-1/2) + (decf e) + (setf x (- (scale-float x 1) 1))) + (t + (decf x))) + (setf z (* x x)) + (setf y (* x (/ (* z (poly-eval x p)) + (poly-eval-1 x q)))) + (decf y (scale-float z -1)))) + (values e x y z))))) + (defun dd-%log2 (x) (declare (type double-double-float x) (optimize (speed 3) (space 0) (inhibit-warnings 3))) - (multiple-value-bind (x e) - (decode-float x) - (declare (type double-double-float x) - (type double-float-exponent e)) - (let ((z 0w0) - (y 0w0)) - (declare (type double-double-float z y)) - (cond ((or (> e 2) - (< e -2)) - (cond ((< x sqrt-1/2) - ;; 2*(2*x-1)/(2*x+1) - (decf e) - (setf z (- x 0.5w0)) - (setf y (+ (* 0.5w0 z) 0.5w0))) - (t - ;; 2*(x-1)/(x+1) - (setf z (- x 0.5w0)) - (decf z 0.5w0) - (setf y (+ (* 0.5w0 x) 0.5w0)))) - (setf x (/ z y)) - (setf z (* x x)) - (setf y (* x (/ (* z (poly-eval z r)) - (poly-eval-1 z s))))) - (t - (cond ((< x sqrt-1/2) - (decf e) - (setf x (- (scale-float x 1) 1))) - (t - (decf x))) - (setf z (* x x)) - (setf y (* x (/ (* z (poly-eval x p)) - (poly-eval-1 x q)))) - (decf y (scale-float z -1)))) - ;; Multiply log of fraction by log2(e) and base 2 exponent by 1 - ;; - ;; This sequence of operations is critical - (setf z (* y log2ea)) - (setf z (+ z (* x log2ea))) - (setf z (+ z y)) - (setf z (+ z x)) - (setf z (+ z e)) - z)))) + (multiple-value-bind (e x y z) + (compute-log x) + ;; Multiply log of fraction by log2(e) and base 2 exponent by 1 + ;; + ;; This sequence of operations is critical + (setf z (* y log2ea)) + (incf z (* x log2ea)) + (incf z y) + (incf z x) + (incf z e) + z)) + + (defun dd-%log10 (x) + (declare (type double-double-float x) + (optimize (speed 3) (space 0) + (inhibit-warnings 3))) + (multiple-value-bind (e x y z) + (compute-log x) + ;; Multiply log of fraction by log10(e) and base 2 exponent by log10(2). + ;; + ;; This sequence of operations is critical. + (setf z (* y l10eb)) + (incf z (* x l10eb)) + (incf z (* e l102b)) + (incf z (* y l10ea)) + (incf z (* x l10ea)) + (incf z (* e l102a)) + z)))) ;;; dd-%exp2 ;;; 2^x ----------------------------------------------------------------------- Summary of changes: src/code/irrat-dd.lisp | 135 +++++++++++++++++++++++++++++++++---------------- 1 file changed, 91 insertions(+), 44 deletions(-) hooks/post-receive -- CMU Common Lisp

1 0