<html lang='en'>
<head>
<meta content='text/html; charset=utf-8' http-equiv='Content-Type'>
<title>
GitLab
</title>
</meta>
</head>
<style>
img {
max-width: 100%;
height: auto;
}
p.details {
font-style:italic;
color:#777
}
.footer p {
font-size:small;
color:#777
}
pre.commit-message {
white-space: pre-wrap;
}
.file-stats a {
text-decoration: none;
}
.file-stats .new-file {
color: #090;
}
.file-stats .deleted-file {
color: #B00;
}}
</style>
<body>
<div class='content'>
<h3>Raymond Toy pushed to branch master at <a href="https://gitlab.common-lisp.net/cmucl/cmucl">cmucl / cmucl</a></h3>
<h4>
Commits:
</h4>
<ul>
<li>
<strong><a href="https://gitlab.common-lisp.net/cmucl/cmucl/commit/7aa4fe5779d0282cc32dded5a9039ae0e28d518d">7aa4fe57</a></strong>
<div>
<span>by Raymond Toy</span>
<i>at 2015-06-06T17:26:01Z</i>
</div>
<pre class='commit-message'>Move the double-double functions and transforms to their own file.
compiler/float-tran-dd.lisp:
* Most of the double-double implementation moved here.
compiler/float-tran.lisp:
* Removed most of the double-double implementation.
compiler/loadcom.lisp:
* Load float-tran-dd.
tools/comcom.lisp:
* Compile float-tran-dd.
i18n/local/cmucl.pot:
* Regenerated.</pre>
</li>
<li>
<strong><a href="https://gitlab.common-lisp.net/cmucl/cmucl/commit/a2f0af9423366107068a01df40d3eb59c7b03427">a2f0af94</a></strong>
<div>
<span>by Raymond Toy</span>
<i>at 2015-06-06T17:42:15Z</i>
</div>
<pre class='commit-message'>Move more double-double items from float-tran to float-tran-dd.</pre>
</li>
<li>
<strong><a href="https://gitlab.common-lisp.net/cmucl/cmucl/commit/b6a9e4dd893225c2bf6d09353aa22caa63703d3c">b6a9e4dd</a></strong>
<div>
<span>by Raymond Toy</span>
<i>at 2015-06-06T17:54:55Z</i>
</div>
<pre class='commit-message'>Actually remove the double-double stuff that was moved.
Previously only commented out, so really remove them now.</pre>
</li>
<li>
<strong><a href="https://gitlab.common-lisp.net/cmucl/cmucl/commit/268309aceaab58f9e11bc48da9841735dcefb665">268309ac</a></strong>
<div>
<span>by Raymond Toy</span>
<i>at 2015-06-06T17:55:27Z</i>
</div>
<pre class='commit-message'>Remove #+double-double conditionals.
This file should only be compiled when double-double support is
available.</pre>
</li>
<li>
<strong><a href="https://gitlab.common-lisp.net/cmucl/cmucl/commit/e6da7db142dd0e9ed16cd53c012fd9db7fb15642">e6da7db1</a></strong>
<div>
<span>by Raymond Toy</span>
<i>at 2015-06-06T18:24:23Z</i>
</div>
<pre class='commit-message'>Wrap the entire file in double-double conditional.
This is so that we can always compile and load this file, even if
double-double isn't supported. (But all currently supported
architectures support double-doubles.)</pre>
</li>
<li>
<strong><a href="https://gitlab.common-lisp.net/cmucl/cmucl/commit/fba8332d1273932ab4631d223880b986b2736f7e">fba8332d</a></strong>
<div>
<span>by Raymond Toy</span>
<i>at 2015-06-06T18:24:38Z</i>
</div>
<pre class='commit-message'>Regenerated.</pre>
</li>
<li>
<strong><a href="https://gitlab.common-lisp.net/cmucl/cmucl/commit/b7af30bd57a3ab56e0d1d9c398df3200f9cc45d6">b7af30bd</a></strong>
<div>
<span>by Raymond Toy</span>
<i>at 2015-06-07T03:40:47Z</i>
</div>
<pre class='commit-message'>Merge branch 'rtoy-split-out-dd-math' into 'master'
Split out double-double math routines
Move the double-double transforms and a few other double-double methods from float-tran.lisp to float-tran-dd.lisp
See merge request !1</pre>
</li>
</ul>
<h4>5 changed files:</h4>
<ul>
<li class='file-stats'>
<a href='#diff-0'>
<span class='new-file'>
+
src/compiler/float-tran-dd.lisp
</span>
</a>
</li>
<li class='file-stats'>
<a href='#diff-1'>
src/compiler/float-tran.lisp
</a>
</li>
<li class='file-stats'>
<a href='#diff-2'>
src/compiler/loadcom.lisp
</a>
</li>
<li class='file-stats'>
<a href='#diff-3'>
src/i18n/locale/cmucl.pot
</a>
</li>
<li class='file-stats'>
<a href='#diff-4'>
src/tools/comcom.lisp
</a>
</li>
</ul>
<h4>Changes:</h4>
<li id='diff-0'>
<a href='https://gitlab.common-lisp.net/cmucl/cmucl/compare/65a61bdbf178b151a633825e27761e25fa45442e...b7af30bd57a3ab56e0d1d9c398df3200f9cc45d6#diff-0'>
<strong>
src/compiler/float-tran-dd.lisp
</strong>
</a>
<hr>
<pre class="highlight"><code><span style="color: #000000;background-color: #ffdddd">--- /dev/null
</span><span style="color: #000000;background-color: #ddffdd">+++ b/src/compiler/float-tran-dd.lisp
</span><span style="color: #aaaaaa">@@ -0,0 +1,690 @@
</span><span style="color: #000000;background-color: #ddffdd">+;;; -*- Mode: Lisp; Package: C; Log: code.log -*-
+;;;
+;;; **********************************************************************
+;;; 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: src/compiler/float-tran-dd.lisp $")
+;;;
+;;; **********************************************************************
+;;;
+;;; This file contains floating-point specific transforms for
+;;; double-doubles.
+;;;
+;;; The algorithms contained herein are based on the code written by
+;;; Yozo Hida. See http://www.cs.berkeley.edu/~yozo/ for more
+;;; information.
+
+(in-package "C")
+(intl:textdomain "cmucl")
+
+#+double-double
+(progn
+(defknown %double-double-float (real)
+ double-double-float
+ (movable foldable flushable))
+
+(deftransform float ((n prototype) (* double-double-float) * :when :both)
+ '(%double-double-float n))
+
+(deftransform %double-float ((n) (double-double-float) * :when :both)
+ '(double-double-hi n))
+
+(deftransform %single-float ((n) (double-double-float) * :when :both)
+ '(float (double-double-hi n) 1f0))
+
+(deftransform %double-double-float ((n) (double-double-float) * :when :both)
+ 'n)
+
+#+nil
+(defun %double-double-float (n)
+ (make-double-double-float (float n 1d0) 0d0))
+
+;; Moved to code/float.lisp, because we need this relatively early in
+;; the build process to handle float and real types.
+#+nil
+(defun %double-double-float (n)
+ (typecase n
+ (fixnum
+ (%make-double-double-float (float n 1d0) 0d0))
+ (single-float
+ (%make-double-double-float (float n 1d0) 0d0))
+ (double-float
+ (%make-double-double-float (float n 1d0) 0d0))
+ (double-double-float
+ n)
+ (bignum
+ (bignum:bignum-to-float n 'double-double-float))
+ (ratio
+ (kernel::float-ratio n 'double-double-float))))
+
+(defknown double-double-float-p (t)
+ boolean
+ (movable foldable flushable))
+
+(defknown %make-double-double-float (double-float double-float)
+ double-double-float
+ (movable foldable flushable))
+
+
+(defknown double-double-hi (double-double-float)
+ double-float
+ (movable foldable flushable))
+
+(defknown double-double-lo (double-double-float)
+ double-float
+ (movable foldable flushable))
+
+(deftransform float-sign ((float &optional float2)
+ (double-double-float &optional double-double-float) *)
+ (if float2
+ (let ((temp (gensym)))
+ `(let ((,temp (abs float2)))
+ (if (minusp (float-sign (double-double-hi float)))
+ (- ,temp)
+ ,temp)))
+ '(if (minusp (float-sign (double-double-hi float))) -1w0 1w0)))
+
+(deftransform cis ((x) (double-double-float) *)
+ `(multiple-value-bind (s c)
+ (kernel::dd-%sincos x)
+ (complex c s)))
+
+
+
+(declaim (inline quick-two-sum))
+(defun quick-two-sum (a b)
+ "Computes fl(a+b) and err(a+b), assuming |a| >= |b|"
+ (declare (double-float a b))
+ (let* ((s (+ a b))
+ (e (- b (- s a))))
+ (values s e)))
+
+(declaim (inline two-sum))
+(defun two-sum (a b)
+ "Computes fl(a+b) and err(a+b)"
+ (declare (double-float a b))
+ (let* ((s (+ a b))
+ (v (- s a))
+ (e (+ (- a (- s v))
+ (- b v))))
+ (locally
+ (declare (optimize (inhibit-warnings 3)))
+ (values s e))))
+
+(declaim (maybe-inline add-dd))
+(defun add-dd (a0 a1 b0 b1)
+ "Add the double-double A0,A1 to the double-double B0,B1"
+ (declare (double-float a0 a1 b0 b1)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (multiple-value-bind (s1 s2)
+ (two-sum a0 b0)
+ (declare (double-float s1 s2))
+ (when (float-infinity-p s1)
+ (return-from add-dd (values s1 0d0)))
+ (multiple-value-bind (t1 t2)
+ (two-sum a1 b1)
+ (declare (double-float t1 t2))
+ (incf s2 t1)
+ (multiple-value-bind (s1 s2)
+ (quick-two-sum s1 s2)
+ (declare (double-float s1 s2))
+ (incf s2 t2)
+ (multiple-value-bind (r1 r2)
+ (quick-two-sum s1 s2)
+ (if (and (zerop a0) (zerop b0))
+ ;; Handle sum of signed zeroes here.
+ (values (float-sign (+ a0 b0) 0d0)
+ 0d0)
+ (values r1 r2)))))))
+
+(deftransform + ((a b) (vm::double-double-float vm::double-double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (add-dd (kernel:double-double-hi a) (kernel:double-double-lo a)
+ (kernel:double-double-hi b) (kernel:double-double-lo b))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(declaim (inline quick-two-diff))
+(defun quick-two-diff (a b)
+ "Compute fl(a-b) and err(a-b), assuming |a| >= |b|"
+ (declare (double-float a b))
+ (let ((s (- a b)))
+ (values s (- (- a s) b))))
+
+(declaim (inline two-diff))
+(defun two-diff (a b)
+ "Compute fl(a-b) and err(a-b)"
+ (declare (double-float a b))
+ (let* ((s (- a b))
+ (v (- s a))
+ (e (- (- a (- s v))
+ (+ b v))))
+ (locally
+ (declare (optimize (inhibit-warnings 3)))
+ (values s e))))
+
+(declaim (maybe-inline sub-dd))
+(defun sub-dd (a0 a1 b0 b1)
+ "Subtract the double-double B0,B1 from A0,A1"
+ (declare (double-float a0 a1 b0 b1)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (multiple-value-bind (s1 s2)
+ (two-diff a0 b0)
+ (declare (double-float s2))
+ (when (float-infinity-p s1)
+ (return-from sub-dd (values s1 0d0)))
+ (multiple-value-bind (t1 t2)
+ (two-diff a1 b1)
+ (incf s2 t1)
+ (multiple-value-bind (s1 s2)
+ (quick-two-sum s1 s2)
+ (declare (double-float s2))
+ (incf s2 t2)
+ (multiple-value-bind (r1 r2)
+ (quick-two-sum s1 s2)
+ (if (and (zerop a0) (zerop b0))
+ (values (float-sign (- a0 b0) 0d0)
+ 0d0)
+ (values r1 r2)))))))
+
+(declaim (maybe-inline sub-d-dd))
+(defun sub-d-dd (a b0 b1)
+ "Compute double-double = double - double-double"
+ (declare (double-float a b0 b1)
+ (optimize (speed 3) (safety 0)
+ (inhibit-warnings 3)))
+ (multiple-value-bind (s1 s2)
+ (two-diff a b0)
+ (declare (double-float s2))
+ (when (float-infinity-p s1)
+ (return-from sub-d-dd (values s1 0d0)))
+ (decf s2 b1)
+ (multiple-value-bind (r1 r2)
+ (quick-two-sum s1 s2)
+ (if (and (zerop a) (zerop b0))
+ (values (float-sign (- a b0) 0d0) 0d0)
+ (values r1 r2)))))
+
+(declaim (maybe-inline sub-dd-d))
+(defun sub-dd-d (a0 a1 b)
+ "Subtract the double B from the double-double A0,A1"
+ (declare (double-float a0 a1 b)
+ (optimize (speed 3) (safety 0)
+ (inhibit-warnings 3)))
+ (multiple-value-bind (s1 s2)
+ (two-diff a0 b)
+ (declare (double-float s2))
+ (when (float-infinity-p s1)
+ (return-from sub-dd-d (values s1 0d0)))
+ (incf s2 a1)
+ (multiple-value-bind (r1 r2)
+ (quick-two-sum s1 s2)
+ (if (and (zerop a0) (zerop b))
+ (values (float-sign (- a0 b) 0d0) 0d0)
+ (values r1 r2)))))
+
+(deftransform - ((a b) (vm::double-double-float vm::double-double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (sub-dd (kernel:double-double-hi a) (kernel:double-double-lo a)
+ (kernel:double-double-hi b) (kernel:double-double-lo b))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(deftransform - ((a b) (double-float vm::double-double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (sub-d-dd a
+ (kernel:double-double-hi b) (kernel:double-double-lo b))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(deftransform - ((a b) (vm::double-double-float double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (sub-dd-d (kernel:double-double-hi a) (kernel:double-double-lo a)
+ b)
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(declaim (maybe-inline split))
+;; See Listing 2.6: Mul12 in "CR-LIBM: A library of correctly rounded
+;; elementary functions in double-precision". Also known as Dekker's
+;; algorithm.
+(defun split (a)
+ "Split the double-float number a into a-hi and a-lo such that a =
+ a-hi + a-lo and a-hi contains the upper 26 significant bits of a and
+ a-lo contains the lower 26 bits."
+ (declare (double-float a))
+ (let* ((tmp (* a (+ 1 (expt 2 27))))
+ (a-hi (- tmp (- tmp a)))
+ (a-lo (- a a-hi)))
+ (values a-hi a-lo)))
+
+;; Values used for scaling in two-prod. These are used to determine
+;; if SPLIT might overflow so the value (and result) can be scaled to
+;; prevent overflow.
+(defconstant +two970+
+ (scale-float 1d0 970))
+
+(defconstant +two53+
+ (scale-float 1d0 53))
+
+(defconstant +two-53+
+ (scale-float 1d0 -53))
+
+(declaim (inline two-prod))
+
+;; This is essentially the algorithm given by Listing 2.7 Mul12Cond
+;; given in "CR-LIBM: A library of correctly rounded elementary
+;; functions in double-precision".
+#-ppc
+(defun two-prod (a b)
+ _N"Compute fl(a*b) and err(a*b)"
+ (declare (double-float a b)
+ (optimize (speed 3)))
+ ;; If the numbers are too big, scale them done so SPLIT doesn't overflow.
+ (multiple-value-bind (aa bb)
+ (values (if (> a +two970+)
+ (* a +two-53+)
+ a)
+ (if (> b +two970+)
+ (* b +two-53+)
+ b))
+ (let ((p (* aa bb)))
+ (declare (double-float p)
+ (inline split))
+ (multiple-value-bind (aa-hi aa-lo)
+ (split aa)
+ ;;(format t "aa-hi, aa-lo = ~S ~S~%" aa-hi aa-lo)
+ (multiple-value-bind (bb-hi bb-lo)
+ (split bb)
+ ;;(format t "bb-hi, bb-lo = ~S ~S~%" bb-hi bb-lo)
+ (let ((e (+ (+ (- (* aa-hi bb-hi) p)
+ (* aa-hi bb-lo)
+ (* aa-lo bb-hi))
+ (* aa-lo bb-lo))))
+ (declare (double-float e))
+ (locally
+ (declare (optimize (inhibit-warnings 3)))
+ ;; If the numbers was scaled down, we need to scale the
+ ;; result back up.
+ (when (> a +two970+)
+ (setf p (* p +two53+)
+ e (* e +two53+)))
+ (when (> b +two970+)
+ (setf p (* p +two53+)
+ e (* e +two53+)))
+ (values p e))))))))
+
+#+ppc
+(defun two-prod (a b)
+ _N"Compute fl(a*b) and err(a*b)"
+ (declare (double-float a b))
+ ;; PPC has a fused multiply-subtract instruction that can be used
+ ;; here, so use it.
+ (let* ((p (* a b))
+ (err (vm::fused-multiply-subtract a b p)))
+ (values p err)))
+
+(declaim (inline two-sqr))
+#-ppc
+(defun two-sqr (a)
+ _N"Compute fl(a*a) and err(a*b). This is a more efficient
+ implementation of two-prod"
+ (declare (double-float a))
+ (let ((q (* a a)))
+ (multiple-value-bind (a-hi a-lo)
+ (split a)
+ (locally
+ (declare (optimize (inhibit-warnings 3)))
+ (values q (+ (+ (- (* a-hi a-hi) q)
+ (* 2 a-hi a-lo))
+ (* a-lo a-lo)))))))
+(defun two-sqr (a)
+ _N"Compute fl(a*a) and err(a*b). This is a more efficient
+ implementation of two-prod"
+ (declare (double-float a))
+ (let ((aa (if (> a +two970+)
+ (* a +two-53+)
+ a)))
+ (let ((q (* aa aa)))
+ (declare (double-float q)
+ (inline split))
+ (multiple-value-bind (a-hi a-lo)
+ (split aa)
+ (locally
+ (declare (optimize (inhibit-warnings 3)))
+ (let ((e (+ (+ (- (* a-hi a-hi) q)
+ (* 2 a-hi a-lo))
+ (* a-lo a-lo))))
+ (if (> a +two970+)
+ (values (* q +two53+)
+ (* e +two53+))
+ (values q e))))))))
+
+#+ppc
+(defun two-sqr (a)
+ _N"Compute fl(a*a) and err(a*b). This is a more efficient
+ implementation of two-prod"
+ (declare (double-float a))
+ (let ((q (* a a)))
+ (values q (vm::fused-multiply-subtract a a q))))
+
+(declaim (maybe-inline mul-dd-d))
+(defun mul-dd-d (a0 a1 b)
+ (declare (double-float a0 a1 b)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (multiple-value-bind (p1 p2)
+ (two-prod a0 b)
+ (declare (double-float p2))
+ (when (float-infinity-p p1)
+ (return-from mul-dd-d (values p1 0d0)))
+ ;;(format t "mul-dd-d p1,p2 = ~A ~A~%" p1 p2)
+ (incf p2 (* a1 b))
+ ;;(format t "mul-dd-d p2 = ~A~%" p2)
+ (multiple-value-bind (r1 r2)
+ (quick-two-sum p1 p2)
+ (when (zerop r1)
+ (setf r1 (float-sign p1 0d0))
+ (setf r2 p1))
+ (values r1 r2))))
+
+(declaim (maybe-inline mul-dd))
+(defun mul-dd (a0 a1 b0 b1)
+ "Multiply the double-double A0,A1 with B0,B1"
+ (declare (double-float a0 a1 b0 b1)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (multiple-value-bind (p1 p2)
+ (two-prod a0 b0)
+ (declare (double-float p1 p2))
+ (when (float-infinity-p p1)
+ (return-from mul-dd (values p1 0d0)))
+ (incf p2 (* a0 b1))
+ (incf p2 (* a1 b0))
+ (multiple-value-bind (r1 r2)
+ (quick-two-sum p1 p2)
+ (if (zerop r1)
+ (values (float-sign p1 0d0) 0d0)
+ (values r1 r2)))))
+
+(declaim (maybe-inline add-dd-d))
+(defun add-dd-d (a0 a1 b)
+ "Add the double-double A0,A1 to the double B"
+ (declare (double-float a0 a1 b)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (multiple-value-bind (s1 s2)
+ (two-sum a0 b)
+ (declare (double-float s1 s2))
+ (when (float-infinity-p s1)
+ (return-from add-dd-d (values s1 0d0)))
+ (incf s2 a1)
+ (multiple-value-bind (r1 r2)
+ (quick-two-sum s1 s2)
+ (if (and (zerop a0) (zerop b))
+ (values (float-sign (+ a0 b) 0d0) 0d0)
+ (values r1 r2)))))
+
+(declaim (maybe-inline sqr-dd))
+(defun sqr-dd (a0 a1)
+ (declare (double-float a0 a1)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (multiple-value-bind (p1 p2)
+ (two-sqr a0)
+ (declare (double-float p1 p2))
+ (incf p2 (* 2 a0 a1))
+ ;; Hida's version of sqr (qd-2.1.210) has the following line for
+ ;; the sqr function. But if you compare this with mul-dd, this
+ ;; doesn't exist there, and if you leave it in, it produces
+ ;; results that are different from using mul-dd to square a value.
+ #+nil
+ (incf p2 (* a1 a1))
+ (quick-two-sum p1 p2)))
+
+(deftransform + ((a b) (vm::double-double-float (or integer single-float double-float))
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (add-dd-d (kernel:double-double-hi a) (kernel:double-double-lo a)
+ (float b 1d0))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(deftransform + ((a b) ((or integer single-float double-float) vm::double-double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (add-dd-d (kernel:double-double-hi b) (kernel:double-double-lo b)
+ (float a 1d0))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(deftransform * ((a b) (vm::double-double-float vm::double-double-float)
+ * :node node)
+ ;; non-const-same-leaf-ref-p is stolen from two-arg-derive-type.
+ (flet ((non-const-same-leaf-ref-p (x y)
+ ;; Just like same-leaf-ref-p, but we don't care if the
+ ;; value of the leaf is constant or not.
+ (declare (type continuation x y))
+ (let ((x-use (continuation-use x))
+ (y-use (continuation-use y)))
+ (and (ref-p x-use)
+ (ref-p y-use)
+ (eq (ref-leaf x-use) (ref-leaf y-use))))))
+ (destructuring-bind (arg1 arg2)
+ (combination-args node)
+ ;; If the two args to * are the same, we square the number
+ ;; instead of multiply. Squaring is simpler than a full
+ ;; multiply.
+ (if (non-const-same-leaf-ref-p arg1 arg2)
+ `(multiple-value-bind (hi lo)
+ (sqr-dd (kernel:double-double-hi a) (kernel:double-double-lo a))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo)))
+ `(multiple-value-bind (hi lo)
+ (mul-dd (kernel:double-double-hi a) (kernel:double-double-lo a)
+ (kernel:double-double-hi b) (kernel:double-double-lo b))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo)))))))
+
+(deftransform * ((a b) (vm::double-double-float (or integer single-float double-float))
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (mul-dd-d (kernel:double-double-hi a) (kernel:double-double-lo a)
+ (float b 1d0))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(deftransform * ((a b) ((or integer single-float double-float) vm::double-double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (mul-dd-d (kernel:double-double-hi b) (kernel:double-double-lo b)
+ (float a 1d0))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(declaim (maybe-inline div-dd))
+(defun div-dd (a0 a1 b0 b1)
+ "Divide the double-double A0,A1 by B0,B1"
+ (declare (double-float a0 a1 b0 b1)
+ (optimize (speed 3)
+ (inhibit-warnings 3))
+ (inline sub-dd))
+ (let ((q1 (/ a0 b0)))
+ (when (float-infinity-p q1)
+ (return-from div-dd (values q1 0d0)))
+ ;; (q1b0, q1b1) = q1*(b0,b1)
+ ;;(format t "q1 = ~A~%" q1)
+ (multiple-value-bind (q1b0 q1b1)
+ (mul-dd-d b0 b1 q1)
+ ;;(format t "q1*b = ~A ~A~%" q1b0 q1b1)
+ (multiple-value-bind (r0 r1)
+ ;; r = a - q1 * b
+ (sub-dd a0 a1 q1b0 q1b1)
+ ;;(format t "r = ~A ~A~%" r0 r1)
+ (let ((q2 (/ r0 b0)))
+ (multiple-value-bind (q2b0 q2b1)
+ (mul-dd-d b0 b1 q2)
+ (multiple-value-bind (r0 r1)
+ ;; r = r - (q2*b)
+ (sub-dd r0 r1 q2b0 q2b1)
+ (declare (ignore r1))
+ (let ((q3 (/ r0 b0)))
+ (multiple-value-bind (q1 q2)
+ (quick-two-sum q1 q2)
+ (add-dd-d q1 q2 q3))))))))))
+
+(declaim (maybe-inline div-dd-d))
+(defun div-dd-d (a0 a1 b)
+ (declare (double-float a0 a1 b)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (let ((q1 (/ a0 b)))
+ ;; q1 = approx quotient
+ ;; Now compute a - q1 * b
+ (multiple-value-bind (p1 p2)
+ (two-prod q1 b)
+ (multiple-value-bind (s e)
+ (two-diff a0 p1)
+ (declare (double-float e))
+ (incf e a1)
+ (decf e p2)
+ ;; Next approx
+ (let ((q2 (/ (+ s e) b)))
+ (quick-two-sum q1 q2))))))
+
+(deftransform / ((a b) (vm::double-double-float vm::double-double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (div-dd (kernel:double-double-hi a) (kernel:double-double-lo a)
+ (kernel:double-double-hi b) (kernel:double-double-lo b))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(deftransform / ((a b) (vm::double-double-float (or integer single-float double-float))
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (div-dd-d (kernel:double-double-hi a) (kernel:double-double-lo a)
+ (float b 1d0))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(declaim (inline sqr-d))
+(defun sqr-d (a)
+ "Square"
+ (declare (double-float a)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (two-sqr a))
+
+(declaim (inline mul-d-d))
+(defun mul-d-d (a b)
+ (two-prod a b))
+
+(declaim (maybe-inline sqrt-dd))
+(defun sqrt-dd (a0 a1)
+ (declare (type (double-float 0d0) a0)
+ (double-float a1)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ ;; Strategy: Use Karp's trick: if x is an approximation to sqrt(a),
+ ;; then
+ ;;
+ ;; y = a*x + (a-(a*x)^2)*x/2
+ ;;
+ ;; is an approximation that is accurate to twice the accuracy of x.
+ ;; Also, the multiplication (a*x) and [-]*x can be done with only
+ ;; half the precision.
+ (if (and (zerop a0) (zerop a1))
+ (values a0 a1)
+ (let* ((x (/ (sqrt a0)))
+ (ax (* a0 x)))
+ (multiple-value-bind (s0 s1)
+ (sqr-d ax)
+ (multiple-value-bind (s2)
+ (sub-dd a0 a1 s0 s1)
+ (multiple-value-bind (p0 p1)
+ (mul-d-d s2 (* x 0.5d0))
+ (add-dd-d p0 p1 ax)))))))
+
+(deftransform sqrt ((a) ((vm::double-double-float 0w0))
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (sqrt-dd (kernel:double-double-hi a) (kernel:double-double-lo a))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(declaim (inline neg-dd))
+(defun neg-dd (a0 a1)
+ (declare (double-float a0 a1)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (values (- a0) (- a1)))
+
+(declaim (inline abs-dd))
+(defun abs-dd (a0 a1)
+ (declare (double-float a0 a1)
+ (optimize (speed 3)
+ (inhibit-warnings 3)))
+ (if (minusp a0)
+ (neg-dd a0 a1)
+ (values a0 a1)))
+
+(deftransform abs ((a) (vm::double-double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (abs-dd (kernel:double-double-hi a) (kernel:double-double-lo a))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(deftransform %negate ((a) (vm::double-double-float)
+ * :node node)
+ `(multiple-value-bind (hi lo)
+ (neg-dd (kernel:double-double-hi a) (kernel:double-double-lo a))
+ (truly-the ,(type-specifier (node-derived-type node))
+ (kernel:%make-double-double-float hi lo))))
+
+(declaim (inline dd=))
+(defun dd= (a0 a1 b0 b1)
+ (and (= a0 b0)
+ (= a1 b1)))
+
+(declaim (inline dd<))
+(defun dd< (a0 a1 b0 b1)
+ (or (< a0 b0)
+ (and (= a0 b0)
+ (< a1 b1))))
+
+(declaim (inline dd>))
+(defun dd> (a0 a1 b0 b1)
+ (or (> a0 b0)
+ (and (= a0 b0)
+ (> a1 b1))))
+
+(deftransform = ((a b) (vm::double-double-float vm::double-double-float) *)
+ `(dd= (kernel:double-double-hi a)
+ (kernel:double-double-lo a)
+ (kernel:double-double-hi b)
+ (kernel:double-double-lo b)))
+
+
+(deftransform < ((a b) (vm::double-double-float vm::double-double-float) *)
+ `(dd< (kernel:double-double-hi a)
+ (kernel:double-double-lo a)
+ (kernel:double-double-hi b)
+ (kernel:double-double-lo b)))
+
+
+(deftransform > ((a b) (vm::double-double-float vm::double-double-float) *)
+ `(dd> (kernel:double-double-hi a)
+ (kernel:double-double-lo a)
+ (kernel:double-double-hi b)
+ (kernel:double-double-lo b)))
+) ; end progn
</span></code></pre>
<br>
</li>
<li id='diff-1'>
<a href='https://gitlab.common-lisp.net/cmucl/cmucl/compare/65a61bdbf178b151a633825e27761e25fa45442e...b7af30bd57a3ab56e0d1d9c398df3200f9cc45d6#diff-1'>
<strong>
src/compiler/float-tran.lisp
</strong>
</a>
<hr>
<pre class="highlight"><code><span style="color: #000000;background-color: #ffdddd">--- a/src/compiler/float-tran.lisp
</span><span style="color: #000000;background-color: #ddffdd">+++ b/src/compiler/float-tran.lisp
</span><span style="color: #aaaaaa">@@ -38,47 +38,6 @@
</span> (deftransform %double-float ((n) (double-float) * :when :both)
'n)
-#+double-double
-(progn
-(defknown %double-double-float (real)
<span style="color: #000000;background-color: #ffdddd">- double-double-float
- (movable foldable flushable))
</span>-
-(deftransform float ((n prototype) (* double-double-float) * :when :both)
<span style="color: #000000;background-color: #ffdddd">- '(%double-double-float n))
</span>-
-(deftransform %double-float ((n) (double-double-float) * :when :both)
<span style="color: #000000;background-color: #ffdddd">- '(double-double-hi n))
</span>-
-(deftransform %single-float ((n) (double-double-float) * :when :both)
<span style="color: #000000;background-color: #ffdddd">- '(float (double-double-hi n) 1f0))
</span>-
-(deftransform %double-double-float ((n) (double-double-float) * :when :both)
<span style="color: #000000;background-color: #ffdddd">- 'n)
</span>-
-#+nil
-(defun %double-double-float (n)
<span style="color: #000000;background-color: #ffdddd">- (make-double-double-float (float n 1d0) 0d0))
</span>-
-;; Moved to code/float.lisp, because we need this relatively early in
-;; the build process to handle float and real types.
-#+nil
-(defun %double-double-float (n)
<span style="color: #000000;background-color: #ffdddd">- (typecase n
- (fixnum
- (%make-double-double-float (float n 1d0) 0d0))
- (single-float
- (%make-double-double-float (float n 1d0) 0d0))
- (double-float
- (%make-double-double-float (float n 1d0) 0d0))
- (double-double-float
- n)
- (bignum
- (bignum:bignum-to-float n 'double-double-float))
- (ratio
- (kernel::float-ratio n 'double-double-float))))
</span>-); progn
-
(defknown %complex-single-float (number) (complex single-float)
(movable foldable flushable))
(defknown %complex-double-float (number) (complex double-float)
<span style="color: #aaaaaa">@@ -364,27 +323,6 @@
</span> (values (signed-byte 32) (unsigned-byte 32))
(movable foldable flushable))
-#+double-double
-(progn
-(defknown double-double-float-p (t)
<span style="color: #000000;background-color: #ffdddd">- boolean
- (movable foldable flushable))
</span>-
-(defknown %make-double-double-float (double-float double-float)
<span style="color: #000000;background-color: #ffdddd">- double-double-float
- (movable foldable flushable))
</span>-
-
-(defknown double-double-hi (double-double-float)
<span style="color: #000000;background-color: #ffdddd">- double-float
- (movable foldable flushable))
</span>-
-(defknown double-double-lo (double-double-float)
<span style="color: #000000;background-color: #ffdddd">- double-float
- (movable foldable flushable))
</span>-
-) ; progn
-
(deftransform float-sign ((float &optional float2)
(single-float &optional single-float) *)
(if float2
<span style="color: #aaaaaa">@@ -401,18 +339,6 @@
</span> (if (minusp (double-float-high-bits float)) (- ,temp) ,temp)))
'(if (minusp (double-float-high-bits float)) -1d0 1d0)))
-#+double-double
-(deftransform float-sign ((float &optional float2)
- (double-double-float &optional double-double-float) *)
<span style="color: #000000;background-color: #ffdddd">- (if float2
- (let ((temp (gensym)))
</span>- `(let ((,temp (abs float2)))
- (if (minusp (float-sign (double-double-hi float)))
- (- ,temp)
- ,temp)))
<span style="color: #000000;background-color: #ffdddd">- '(if (minusp (float-sign (double-double-hi float))) -1w0 1w0)))
</span>-
<span style="color: #000000;background-color: #ffdddd">-
</span>
;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, SCALE-FLOAT:
;;;
<span style="color: #aaaaaa">@@ -778,13 +704,6 @@
</span> (%sincos x)
(complex c s)))
-#+double-double
-(deftransform cis ((x) (double-double-float) *)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (s c)
- (kernel::dd-%sincos x)
- (complex c s)))
</span>-
-
;;; The argument range is limited on the x86 FP trig. functions. A
;;; post-test can detect a failure (and load a suitable result), but
;;; this test is avoided if possible.
<span style="color: #aaaaaa">@@ -2170,609 +2089,3 @@
</span> (make-values-type :required (list f
e
s))))
-
-;;; Support for double-double floats
-;;;
-;;; The algorithms contained herein are based on the code written by
-;;; Yozo Hida. See http://www.cs.berkeley.edu/~yozo/ for more
-;;; information.
-
-#+double-double
-(progn
<span style="color: #000000;background-color: #ffdddd">-
</span>-(declaim (inline quick-two-sum))
-(defun quick-two-sum (a b)
<span style="color: #000000;background-color: #ffdddd">- "Computes fl(a+b) and err(a+b), assuming |a| >= |b|"
- (declare (double-float a b))
- (let* ((s (+ a b))
</span>- (e (- b (- s a))))
<span style="color: #000000;background-color: #ffdddd">- (values s e)))
</span>-
-(declaim (inline two-sum))
-(defun two-sum (a b)
<span style="color: #000000;background-color: #ffdddd">- "Computes fl(a+b) and err(a+b)"
- (declare (double-float a b))
- (let* ((s (+ a b))
</span>- (v (- s a))
- (e (+ (- a (- s v))
- (- b v))))
<span style="color: #000000;background-color: #ffdddd">- (locally
</span>- (declare (optimize (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (values s e))))
</span>-
-(declaim (maybe-inline add-dd))
-(defun add-dd (a0 a1 b0 b1)
<span style="color: #000000;background-color: #ffdddd">- "Add the double-double A0,A1 to the double-double B0,B1"
- (declare (double-float a0 a1 b0 b1)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (s1 s2)
- (two-sum a0 b0)
- (declare (double-float s1 s2))
- (when (float-infinity-p s1)
- (return-from add-dd (values s1 0d0)))
- (multiple-value-bind (t1 t2)
</span>- (two-sum a1 b1)
<span style="color: #000000;background-color: #ffdddd">- (declare (double-float t1 t2))
- (incf s2 t1)
- (multiple-value-bind (s1 s2)
</span>- (quick-two-sum s1 s2)
- (declare (double-float s1 s2))
- (incf s2 t2)
- (multiple-value-bind (r1 r2)
- (quick-two-sum s1 s2)
- (if (and (zerop a0) (zerop b0))
- ;; Handle sum of signed zeroes here.
- (values (float-sign (+ a0 b0) 0d0)
- 0d0)
- (values r1 r2)))))))
-
-(deftransform + ((a b) (vm::double-double-float vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (add-dd (kernel:double-double-hi a) (kernel:double-double-lo a)
</span>- (kernel:double-double-hi b) (kernel:double-double-lo b))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(declaim (inline quick-two-diff))
-(defun quick-two-diff (a b)
<span style="color: #000000;background-color: #ffdddd">- "Compute fl(a-b) and err(a-b), assuming |a| >= |b|"
- (declare (double-float a b))
- (let ((s (- a b)))
- (values s (- (- a s) b))))
</span>-
-(declaim (inline two-diff))
-(defun two-diff (a b)
<span style="color: #000000;background-color: #ffdddd">- "Compute fl(a-b) and err(a-b)"
- (declare (double-float a b))
- (let* ((s (- a b))
</span>- (v (- s a))
- (e (- (- a (- s v))
- (+ b v))))
<span style="color: #000000;background-color: #ffdddd">- (locally
</span>- (declare (optimize (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (values s e))))
</span>-
-(declaim (maybe-inline sub-dd))
-(defun sub-dd (a0 a1 b0 b1)
<span style="color: #000000;background-color: #ffdddd">- "Subtract the double-double B0,B1 from A0,A1"
- (declare (double-float a0 a1 b0 b1)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (s1 s2)
- (two-diff a0 b0)
- (declare (double-float s2))
- (when (float-infinity-p s1)
- (return-from sub-dd (values s1 0d0)))
- (multiple-value-bind (t1 t2)
</span>- (two-diff a1 b1)
<span style="color: #000000;background-color: #ffdddd">- (incf s2 t1)
- (multiple-value-bind (s1 s2)
</span>- (quick-two-sum s1 s2)
- (declare (double-float s2))
- (incf s2 t2)
- (multiple-value-bind (r1 r2)
- (quick-two-sum s1 s2)
- (if (and (zerop a0) (zerop b0))
- (values (float-sign (- a0 b0) 0d0)
- 0d0)
- (values r1 r2)))))))
-
-(declaim (maybe-inline sub-d-dd))
-(defun sub-d-dd (a b0 b1)
<span style="color: #000000;background-color: #ffdddd">- "Compute double-double = double - double-double"
- (declare (double-float a b0 b1)
</span>- (optimize (speed 3) (safety 0)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (s1 s2)
- (two-diff a b0)
- (declare (double-float s2))
- (when (float-infinity-p s1)
- (return-from sub-d-dd (values s1 0d0)))
- (decf s2 b1)
- (multiple-value-bind (r1 r2)
</span>- (quick-two-sum s1 s2)
<span style="color: #000000;background-color: #ffdddd">- (if (and (zerop a) (zerop b0))
</span>- (values (float-sign (- a b0) 0d0) 0d0)
- (values r1 r2)))))
-
-(declaim (maybe-inline sub-dd-d))
-(defun sub-dd-d (a0 a1 b)
<span style="color: #000000;background-color: #ffdddd">- "Subtract the double B from the double-double A0,A1"
- (declare (double-float a0 a1 b)
</span>- (optimize (speed 3) (safety 0)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (s1 s2)
- (two-diff a0 b)
- (declare (double-float s2))
- (when (float-infinity-p s1)
- (return-from sub-dd-d (values s1 0d0)))
- (incf s2 a1)
- (multiple-value-bind (r1 r2)
</span>- (quick-two-sum s1 s2)
<span style="color: #000000;background-color: #ffdddd">- (if (and (zerop a0) (zerop b))
</span>- (values (float-sign (- a0 b) 0d0) 0d0)
- (values r1 r2)))))
-
-(deftransform - ((a b) (vm::double-double-float vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (sub-dd (kernel:double-double-hi a) (kernel:double-double-lo a)
</span>- (kernel:double-double-hi b) (kernel:double-double-lo b))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(deftransform - ((a b) (double-float vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (sub-d-dd a
</span>- (kernel:double-double-hi b) (kernel:double-double-lo b))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(deftransform - ((a b) (vm::double-double-float double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (sub-dd-d (kernel:double-double-hi a) (kernel:double-double-lo a)
</span>- b)
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(declaim (maybe-inline split))
-;; See Listing 2.6: Mul12 in "CR-LIBM: A library of correctly rounded
-;; elementary functions in double-precision". Also known as Dekker's
-;; algorithm.
-(defun split (a)
<span style="color: #000000;background-color: #ffdddd">- "Split the double-float number a into a-hi and a-lo such that a =
- a-hi + a-lo and a-hi contains the upper 26 significant bits of a and
- a-lo contains the lower 26 bits."
- (declare (double-float a))
- (let* ((tmp (* a (+ 1 (expt 2 27))))
</span>- (a-hi (- tmp (- tmp a)))
- (a-lo (- a a-hi)))
<span style="color: #000000;background-color: #ffdddd">- (values a-hi a-lo)))
</span>-
-;; Values used for scaling in two-prod. These are used to determine
-;; if SPLIT might overflow so the value (and result) can be scaled to
-;; prevent overflow.
-(defconstant +two970+
<span style="color: #000000;background-color: #ffdddd">- (scale-float 1d0 970))
</span>-
-(defconstant +two53+
<span style="color: #000000;background-color: #ffdddd">- (scale-float 1d0 53))
</span>-
-(defconstant +two-53+
<span style="color: #000000;background-color: #ffdddd">- (scale-float 1d0 -53))
</span>-
-(declaim (inline two-prod))
-
-;; This is essentially the algorithm given by Listing 2.7 Mul12Cond
-;; given in "CR-LIBM: A library of correctly rounded elementary
-;; functions in double-precision".
-#-ppc
-(defun two-prod (a b)
<span style="color: #000000;background-color: #ffdddd">- _N"Compute fl(a*b) and err(a*b)"
- (declare (double-float a b)
</span>- (optimize (speed 3)))
<span style="color: #000000;background-color: #ffdddd">- ;; If the numbers are too big, scale them done so SPLIT doesn't overflow.
- (multiple-value-bind (aa bb)
- (values (if (> a +two970+)
</span>- (* a +two-53+)
- a)
- (if (> b +two970+)
- (* b +two-53+)
- b))
<span style="color: #000000;background-color: #ffdddd">- (let ((p (* aa bb)))
- (declare (double-float p)
</span>- (inline split))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (aa-hi aa-lo)
</span>- (split aa)
- ;;(format t "aa-hi, aa-lo = ~S ~S~%" aa-hi aa-lo)
- (multiple-value-bind (bb-hi bb-lo)
- (split bb)
- ;;(format t "bb-hi, bb-lo = ~S ~S~%" bb-hi bb-lo)
- (let ((e (+ (+ (- (* aa-hi bb-hi) p)
- (* aa-hi bb-lo)
- (* aa-lo bb-hi))
- (* aa-lo bb-lo))))
- (declare (double-float e))
- (locally
- (declare (optimize (inhibit-warnings 3)))
- ;; If the numbers was scaled down, we need to scale the
- ;; result back up.
- (when (> a +two970+)
- (setf p (* p +two53+)
- e (* e +two53+)))
- (when (> b +two970+)
- (setf p (* p +two53+)
- e (* e +two53+)))
- (values p e))))))))
-
-#+ppc
-(defun two-prod (a b)
<span style="color: #000000;background-color: #ffdddd">- _N"Compute fl(a*b) and err(a*b)"
- (declare (double-float a b))
- ;; PPC has a fused multiply-subtract instruction that can be used
- ;; here, so use it.
- (let* ((p (* a b))
</span>- (err (vm::fused-multiply-subtract a b p)))
<span style="color: #000000;background-color: #ffdddd">- (values p err)))
</span>-
-(declaim (inline two-sqr))
-#-ppc
-(defun two-sqr (a)
<span style="color: #000000;background-color: #ffdddd">- _N"Compute fl(a*a) and err(a*b). This is a more efficient
- implementation of two-prod"
- (declare (double-float a))
- (let ((q (* a a)))
- (multiple-value-bind (a-hi a-lo)
</span>- (split a)
<span style="color: #000000;background-color: #ffdddd">- (locally
</span>- (declare (optimize (inhibit-warnings 3)))
- (values q (+ (+ (- (* a-hi a-hi) q)
- (* 2 a-hi a-lo))
- (* a-lo a-lo)))))))
-(defun two-sqr (a)
<span style="color: #000000;background-color: #ffdddd">- _N"Compute fl(a*a) and err(a*b). This is a more efficient
- implementation of two-prod"
- (declare (double-float a))
- (let ((aa (if (> a +two970+)
</span>- (* a +two-53+)
- a)))
<span style="color: #000000;background-color: #ffdddd">- (let ((q (* aa aa)))
- (declare (double-float q)
</span>- (inline split))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (a-hi a-lo)
</span>- (split aa)
- (locally
- (declare (optimize (inhibit-warnings 3)))
- (let ((e (+ (+ (- (* a-hi a-hi) q)
- (* 2 a-hi a-lo))
- (* a-lo a-lo))))
- (if (> a +two970+)
- (values (* q +two53+)
- (* e +two53+))
- (values q e))))))))
-
-#+ppc
-(defun two-sqr (a)
<span style="color: #000000;background-color: #ffdddd">- _N"Compute fl(a*a) and err(a*b). This is a more efficient
- implementation of two-prod"
- (declare (double-float a))
- (let ((q (* a a)))
- (values q (vm::fused-multiply-subtract a a q))))
</span>-
-(declaim (maybe-inline mul-dd-d))
-(defun mul-dd-d (a0 a1 b)
<span style="color: #000000;background-color: #ffdddd">- (declare (double-float a0 a1 b)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (p1 p2)
- (two-prod a0 b)
- (declare (double-float p2))
- (when (float-infinity-p p1)
- (return-from mul-dd-d (values p1 0d0)))
- ;;(format t "mul-dd-d p1,p2 = ~A ~A~%" p1 p2)
- (incf p2 (* a1 b))
- ;;(format t "mul-dd-d p2 = ~A~%" p2)
- (multiple-value-bind (r1 r2)
</span>- (quick-two-sum p1 p2)
<span style="color: #000000;background-color: #ffdddd">- (when (zerop r1)
</span>- (setf r1 (float-sign p1 0d0))
- (setf r2 p1))
<span style="color: #000000;background-color: #ffdddd">- (values r1 r2))))
</span>-
-(declaim (maybe-inline mul-dd))
-(defun mul-dd (a0 a1 b0 b1)
<span style="color: #000000;background-color: #ffdddd">- "Multiply the double-double A0,A1 with B0,B1"
- (declare (double-float a0 a1 b0 b1)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (p1 p2)
- (two-prod a0 b0)
- (declare (double-float p1 p2))
- (when (float-infinity-p p1)
- (return-from mul-dd (values p1 0d0)))
- (incf p2 (* a0 b1))
- (incf p2 (* a1 b0))
- (multiple-value-bind (r1 r2)
</span>- (quick-two-sum p1 p2)
<span style="color: #000000;background-color: #ffdddd">- (if (zerop r1)
</span>- (values (float-sign p1 0d0) 0d0)
- (values r1 r2)))))
-
-(declaim (maybe-inline add-dd-d))
-(defun add-dd-d (a0 a1 b)
<span style="color: #000000;background-color: #ffdddd">- "Add the double-double A0,A1 to the double B"
- (declare (double-float a0 a1 b)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (s1 s2)
- (two-sum a0 b)
- (declare (double-float s1 s2))
- (when (float-infinity-p s1)
- (return-from add-dd-d (values s1 0d0)))
- (incf s2 a1)
- (multiple-value-bind (r1 r2)
</span>- (quick-two-sum s1 s2)
<span style="color: #000000;background-color: #ffdddd">- (if (and (zerop a0) (zerop b))
</span>- (values (float-sign (+ a0 b) 0d0) 0d0)
- (values r1 r2)))))
-
-(declaim (maybe-inline sqr-dd))
-(defun sqr-dd (a0 a1)
<span style="color: #000000;background-color: #ffdddd">- (declare (double-float a0 a1)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (p1 p2)
- (two-sqr a0)
- (declare (double-float p1 p2))
- (incf p2 (* 2 a0 a1))
- ;; Hida's version of sqr (qd-2.1.210) has the following line for
- ;; the sqr function. But if you compare this with mul-dd, this
- ;; doesn't exist there, and if you leave it in, it produces
- ;; results that are different from using mul-dd to square a value.
- #+nil
- (incf p2 (* a1 a1))
- (quick-two-sum p1 p2)))
</span>-
-(deftransform + ((a b) (vm::double-double-float (or integer single-float double-float))
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (add-dd-d (kernel:double-double-hi a) (kernel:double-double-lo a)
</span>- (float b 1d0))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(deftransform + ((a b) ((or integer single-float double-float) vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (add-dd-d (kernel:double-double-hi b) (kernel:double-double-lo b)
</span>- (float a 1d0))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(deftransform * ((a b) (vm::double-double-float vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- ;; non-const-same-leaf-ref-p is stolen from two-arg-derive-type.
- (flet ((non-const-same-leaf-ref-p (x y)
</span>- ;; Just like same-leaf-ref-p, but we don't care if the
- ;; value of the leaf is constant or not.
- (declare (type continuation x y))
- (let ((x-use (continuation-use x))
- (y-use (continuation-use y)))
- (and (ref-p x-use)
- (ref-p y-use)
- (eq (ref-leaf x-use) (ref-leaf y-use))))))
<span style="color: #000000;background-color: #ffdddd">- (destructuring-bind (arg1 arg2)
</span>- (combination-args node)
<span style="color: #000000;background-color: #ffdddd">- ;; If the two args to * are the same, we square the number
- ;; instead of multiply. Squaring is simpler than a full
- ;; multiply.
- (if (non-const-same-leaf-ref-p arg1 arg2)
</span>- `(multiple-value-bind (hi lo)
- (sqr-dd (kernel:double-double-hi a) (kernel:double-double-lo a))
- (truly-the ,(type-specifier (node-derived-type node))
- (kernel:%make-double-double-float hi lo)))
- `(multiple-value-bind (hi lo)
- (mul-dd (kernel:double-double-hi a) (kernel:double-double-lo a)
- (kernel:double-double-hi b) (kernel:double-double-lo b))
- (truly-the ,(type-specifier (node-derived-type node))
- (kernel:%make-double-double-float hi lo)))))))
-
-(deftransform * ((a b) (vm::double-double-float (or integer single-float double-float))
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (mul-dd-d (kernel:double-double-hi a) (kernel:double-double-lo a)
</span>- (float b 1d0))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(deftransform * ((a b) ((or integer single-float double-float) vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (mul-dd-d (kernel:double-double-hi b) (kernel:double-double-lo b)
</span>- (float a 1d0))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(declaim (maybe-inline div-dd))
-(defun div-dd (a0 a1 b0 b1)
<span style="color: #000000;background-color: #ffdddd">- "Divide the double-double A0,A1 by B0,B1"
- (declare (double-float a0 a1 b0 b1)
</span>- (optimize (speed 3)
- (inhibit-warnings 3))
- (inline sub-dd))
<span style="color: #000000;background-color: #ffdddd">- (let ((q1 (/ a0 b0)))
- (when (float-infinity-p q1)
- (return-from div-dd (values q1 0d0)))
- ;; (q1b0, q1b1) = q1*(b0,b1)
- ;;(format t "q1 = ~A~%" q1)
- (multiple-value-bind (q1b0 q1b1)
</span>- (mul-dd-d b0 b1 q1)
<span style="color: #000000;background-color: #ffdddd">- ;;(format t "q1*b = ~A ~A~%" q1b0 q1b1)
- (multiple-value-bind (r0 r1)
</span>- ;; r = a - q1 * b
- (sub-dd a0 a1 q1b0 q1b1)
- ;;(format t "r = ~A ~A~%" r0 r1)
- (let ((q2 (/ r0 b0)))
- (multiple-value-bind (q2b0 q2b1)
- (mul-dd-d b0 b1 q2)
- (multiple-value-bind (r0 r1)
- ;; r = r - (q2*b)
- (sub-dd r0 r1 q2b0 q2b1)
- (declare (ignore r1))
- (let ((q3 (/ r0 b0)))
- (multiple-value-bind (q1 q2)
- (quick-two-sum q1 q2)
- (add-dd-d q1 q2 q3))))))))))
-
-(declaim (maybe-inline div-dd-d))
-(defun div-dd-d (a0 a1 b)
<span style="color: #000000;background-color: #ffdddd">- (declare (double-float a0 a1 b)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (let ((q1 (/ a0 b)))
- ;; q1 = approx quotient
- ;; Now compute a - q1 * b
- (multiple-value-bind (p1 p2)
</span>- (two-prod q1 b)
<span style="color: #000000;background-color: #ffdddd">- (multiple-value-bind (s e)
</span>- (two-diff a0 p1)
- (declare (double-float e))
- (incf e a1)
- (decf e p2)
- ;; Next approx
- (let ((q2 (/ (+ s e) b)))
- (quick-two-sum q1 q2))))))
-
-(deftransform / ((a b) (vm::double-double-float vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (div-dd (kernel:double-double-hi a) (kernel:double-double-lo a)
</span>- (kernel:double-double-hi b) (kernel:double-double-lo b))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(deftransform / ((a b) (vm::double-double-float (or integer single-float double-float))
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (div-dd-d (kernel:double-double-hi a) (kernel:double-double-lo a)
</span>- (float b 1d0))
<span style="color: #000000;background-color: #ffdddd">- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(declaim (inline sqr-d))
-(defun sqr-d (a)
<span style="color: #000000;background-color: #ffdddd">- "Square"
- (declare (double-float a)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (two-sqr a))
</span>-
-(declaim (inline mul-d-d))
-(defun mul-d-d (a b)
<span style="color: #000000;background-color: #ffdddd">- (two-prod a b))
</span>-
-(declaim (maybe-inline sqrt-dd))
-(defun sqrt-dd (a0 a1)
<span style="color: #000000;background-color: #ffdddd">- (declare (type (double-float 0d0) a0)
</span>- (double-float a1)
- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- ;; Strategy: Use Karp's trick: if x is an approximation to sqrt(a),
- ;; then
- ;;
- ;; y = a*x + (a-(a*x)^2)*x/2
- ;;
- ;; is an approximation that is accurate to twice the accuracy of x.
- ;; Also, the multiplication (a*x) and [-]*x can be done with only
- ;; half the precision.
- (if (and (zerop a0) (zerop a1))
- (values a0 a1)
- (let* ((x (/ (sqrt a0)))
</span>- (ax (* a0 x)))
- (multiple-value-bind (s0 s1)
- (sqr-d ax)
- (multiple-value-bind (s2)
- (sub-dd a0 a1 s0 s1)
- (multiple-value-bind (p0 p1)
- (mul-d-d s2 (* x 0.5d0))
- (add-dd-d p0 p1 ax)))))))
-
-(deftransform sqrt ((a) ((vm::double-double-float 0w0))
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (sqrt-dd (kernel:double-double-hi a) (kernel:double-double-lo a))
- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(declaim (inline neg-dd))
-(defun neg-dd (a0 a1)
<span style="color: #000000;background-color: #ffdddd">- (declare (double-float a0 a1)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (values (- a0) (- a1)))
</span>-
-(declaim (inline abs-dd))
-(defun abs-dd (a0 a1)
<span style="color: #000000;background-color: #ffdddd">- (declare (double-float a0 a1)
</span>- (optimize (speed 3)
- (inhibit-warnings 3)))
<span style="color: #000000;background-color: #ffdddd">- (if (minusp a0)
- (neg-dd a0 a1)
- (values a0 a1)))
</span>-
-(deftransform abs ((a) (vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (abs-dd (kernel:double-double-hi a) (kernel:double-double-lo a))
- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(deftransform %negate ((a) (vm::double-double-float)
- * :node node)
<span style="color: #000000;background-color: #ffdddd">- `(multiple-value-bind (hi lo)
- (neg-dd (kernel:double-double-hi a) (kernel:double-double-lo a))
- (truly-the ,(type-specifier (node-derived-type node))
</span>- (kernel:%make-double-double-float hi lo))))
-
-(declaim (inline dd=))
-(defun dd= (a0 a1 b0 b1)
<span style="color: #000000;background-color: #ffdddd">- (and (= a0 b0)
- (= a1 b1)))
-
</span>-(declaim (inline dd<))
-(defun dd< (a0 a1 b0 b1)
<span style="color: #000000;background-color: #ffdddd">- (or (< a0 b0)
- (and (= a0 b0)
</span>- (< a1 b1))))
-
-(declaim (inline dd>))
-(defun dd> (a0 a1 b0 b1)
<span style="color: #000000;background-color: #ffdddd">- (or (> a0 b0)
- (and (= a0 b0)
</span>- (> a1 b1))))
<span style="color: #000000;background-color: #ffdddd">-
</span>-(deftransform = ((a b) (vm::double-double-float vm::double-double-float) *)
<span style="color: #000000;background-color: #ffdddd">- `(dd= (kernel:double-double-hi a)
</span>- (kernel:double-double-lo a)
- (kernel:double-double-hi b)
- (kernel:double-double-lo b)))
-
-
-(deftransform < ((a b) (vm::double-double-float vm::double-double-float) *)
<span style="color: #000000;background-color: #ffdddd">- `(dd< (kernel:double-double-hi a)
</span>- (kernel:double-double-lo a)
- (kernel:double-double-hi b)
- (kernel:double-double-lo b)))
-
-
-(deftransform > ((a b) (vm::double-double-float vm::double-double-float) *)
<span style="color: #000000;background-color: #ffdddd">- `(dd> (kernel:double-double-hi a)
</span>- (kernel:double-double-lo a)
- (kernel:double-double-hi b)
- (kernel:double-double-lo b)))
-
-) ; progn double-double
</code></pre>
<br>
</li>
<li id='diff-2'>
<a href='https://gitlab.common-lisp.net/cmucl/cmucl/compare/65a61bdbf178b151a633825e27761e25fa45442e...b7af30bd57a3ab56e0d1d9c398df3200f9cc45d6#diff-2'>
<strong>
src/compiler/loadcom.lisp
</strong>
</a>
<hr>
<pre class="highlight"><code><span style="color: #000000;background-color: #ffdddd">--- a/src/compiler/loadcom.lisp
</span><span style="color: #000000;background-color: #ddffdd">+++ b/src/compiler/loadcom.lisp
</span><span style="color: #aaaaaa">@@ -32,6 +32,7 @@
</span> (load "vm:vm-typetran")
(load "vm:vm-tran")
(load "c:float-tran")
<span style="color: #000000;background-color: #ddffdd">+(load "c:float-tran-dd")
</span> (load "c:saptran")
(load "c:srctran")
(load "c:locall")
</code></pre>
<br>
</li>
<li id='diff-3'>
<a href='https://gitlab.common-lisp.net/cmucl/cmucl/compare/65a61bdbf178b151a633825e27761e25fa45442e...b7af30bd57a3ab56e0d1d9c398df3200f9cc45d6#diff-3'>
<strong>
src/i18n/locale/cmucl.pot
</strong>
</a>
<hr>
<pre class="highlight"><code><span style="color: #000000;background-color: #ffdddd">--- a/src/i18n/locale/cmucl.pot
</span><span style="color: #000000;background-color: #ddffdd">+++ b/src/i18n/locale/cmucl.pot
</span><span style="color: #aaaaaa">@@ -18776,68 +18776,68 @@ msgstr ""
</span> msgid "Float zero bound ~s not correctly canonicalised?"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Compute fl(a*b) and err(a*b)"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid ""
"Compute fl(a*a) and err(a*b). This is a more efficient\n"
" implementation of two-prod"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Computes fl(a+b) and err(a+b), assuming |a| >= |b|"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Computes fl(a+b) and err(a+b)"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Add the double-double A0,A1 to the double-double B0,B1"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Compute fl(a-b) and err(a-b), assuming |a| >= |b|"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Compute fl(a-b) and err(a-b)"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Subtract the double-double B0,B1 from A0,A1"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Compute double-double = double - double-double"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Subtract the double B from the double-double A0,A1"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid ""
"Split the double-float number a into a-hi and a-lo such that a =\n"
" a-hi + a-lo and a-hi contains the upper 26 significant bits of a and\n"
" a-lo contains the lower 26 bits."
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Multiply the double-double A0,A1 with B0,B1"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Add the double-double A0,A1 to the double B"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Divide the double-double A0,A1 by B0,B1"
msgstr ""
-#: src/compiler/float-tran.lisp
<span style="color: #000000;background-color: #ddffdd">+#: src/compiler/float-tran-dd.lisp
</span> msgid "Square"
msgstr ""
</code></pre>
<br>
</li>
<li id='diff-4'>
<a href='https://gitlab.common-lisp.net/cmucl/cmucl/compare/65a61bdbf178b151a633825e27761e25fa45442e...b7af30bd57a3ab56e0d1d9c398df3200f9cc45d6#diff-4'>
<strong>
src/tools/comcom.lisp
</strong>
</a>
<hr>
<pre class="highlight"><code><span style="color: #000000;background-color: #ffdddd">--- a/src/tools/comcom.lisp
</span><span style="color: #000000;background-color: #ddffdd">+++ b/src/tools/comcom.lisp
</span><span style="color: #aaaaaa">@@ -121,6 +121,7 @@
</span> (comf "target:compiler/typetran" :byte-compile *byte-compile*)
(comf "target:compiler/generic/vm-typetran" :byte-compile *byte-compile*)
(comf "target:compiler/float-tran" :byte-compile *byte-compile*)
<span style="color: #000000;background-color: #ddffdd">+(comf "target:compiler/float-tran-dd" :byte-compile *byte-compile*)
</span> (comf "target:compiler/saptran" :byte-compile *byte-compile*)
(comf "target:compiler/srctran") ;; try
(comf "target:compiler/locall")
</code></pre>
<br>
</li>
</div>
<div class='footer' style='margin-top: 10px;'>
<p>
—
<br>
<a href="https://gitlab.common-lisp.net/cmucl/cmucl/compare/65a61bdbf178b151a633825e27761e25fa45442e...b7af30bd57a3ab56e0d1d9c398df3200f9cc45d6">View it on GitLab</a>
<script type="application/ld+json">{"@context":"http://schema.org","@type":"EmailMessage","action":{"@type":"ViewAction","name":["merge_requests","issues","commit"],"url":"https://gitlab.common-lisp.net/cmucl/cmucl/compare/65a61bdbf178b151a633825e27761e25fa45442e...b7af30bd57a3ab56e0d1d9c398df3200f9cc45d6"}}</script>
</p>
</div>
</body>
</html>