Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: a6ea12d7 by Raymond Toy at 2025-12-13T08:32:38-08:00 Move cdiv-dd-float to float-tran-dd.lisp Move the implementation of cdiv-double-double-float from numbers.lisp to float-tran-dd.lisp. First, we wanted a closure for the constants like we do for cdiv-double-float. However, this requires all the dd routines to be defined and numbers.lisp is compiled too soon. Moving the code to float-tran-dd.lisp makes sure all the functions are ready. Add a deftransform for / to call cdiv-double-double-float. Update two-arg-/ to call cdiv-double-double-float appropriately. - - - - - 8d0c23e1 by Raymond Toy at 2025-12-13T08:35:34-08:00 Import cdiv-dd-float into kernel package and update tests This mean exporting cdiv-double-double-float so that the kernel can import the symbol and use it. Thus, cdiv-dd-float and cdiv-double-float both exist in the kernel package for two-arg-/. Update the float tests to call / instead of cdiv-double-double-float. - - - - - 4 changed files: - src/code/exports.lisp - src/code/numbers.lisp - src/compiler/float-tran-dd.lisp - tests/float.lisp Changes: ===================================== src/code/exports.lisp ===================================== @@ -1669,6 +1669,7 @@ "TRUST-DYNAMIC-EXTENT-DECLARATION-P" "IR2-STACK-ALLOCATE" "%DYNAMIC-EXTENT" "%DYNAMIC-EXTENT-START" "%DYNAMIC-EXTENT-END") + (:export "CDIV-DOUBLE-DOUBLE-FLOAT") ) (defpackage "XREF" (:export "INIT-XREF-DATABASE" @@ -1918,6 +1919,8 @@ (:import-from "LISP" "BOOLEAN") (:import-from "C-CALL" "VOID") (:import-from "CONDITIONS" "PARSE-UNKNOWN-TYPE-SPECIFIER") + #+double-double + (:import-from "C" "CDIV-DOUBLE-DOUBLE-FLOAT") (:shadow "CLASS" "STRUCTURE-CLASS" "BUILT-IN-CLASS" "STANDARD-CLASS" "FIND-CLASS" "CLASS-OF") (:export "%CLASS-LAYOUT" "%CLASS-STATE" "%CLASS-DIRECT-SUPERCLASSES" ===================================== src/code/numbers.lisp ===================================== @@ -752,157 +752,6 @@ (* s (robust-internal x y)))))) -;;; Same algorithm as for doubles, but the constants here are -;;; different. I'm guessing here for the appropriate values for -;;; +eps+, +rmin2+, and +rminscal+. -(defconstant +dd-eps+ (scale-float 1w0 -104)) -(defconstant +dd-rmin+ least-positive-normalized-double-double-float) -(defconstant +dd-rbig+ (/ most-positive-double-double-float 2)) -(defconstant +dd-rmin2+ (scale-float 1w0 -105)) -(defconstant +dd-rminscal+ (scale-float 1w0 102)) -(defconstant +dd-rmax2+ (* +dd-rbig+ +dd-rmin2+)) -(defconstant +dd-be+ (/ 2 (* +dd-eps+ +dd-eps+))) -(defconstant +dd-2/eps+ (/ 2 +dd-eps+)) - -(defun cdiv-double-double-float (x y) - (declare (type (complex vm::double-double-float) x y) - (optimize (speed 3) (safety 0))) - (labels - ((internal-compreal (a b c d r tt) - (declare (vm::double-double-float a b c d r tt)) - ;; Compute the real part of the complex division - ;; (a+ib)/(c+id), assuming |c| <= |d|. r = d/c and tt = 1/(c+d*r). - ;; - ;; The realpart is (a*c+b*d)/(c^2+d^2). - ;; - ;; c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r) - ;; - ;; Then - ;; - ;; (a*c+b*d)/(c^2+d^2) = (a*c+b*d)/(c*(c+d*r)) - ;; = (a + b*d/c)/(c+d*r) - ;; = (a + b*r)/(c + d*r). - ;; - ;; Thus tt = (c + d*r). - (cond ((>= (abs r) +dd-rmin+) - (let ((br (* b r))) - (if (/= br 0) - (/ (+ a br) tt) - ;; b*r underflows. Instead, compute - ;; - ;; (a + b*r)*tt = a*tt + b*tt*r - ;; = a*tt + (b*tt)*r - ;; (a + b*r)/tt = a/tt + b/tt*r - ;; = a*tt + (b*tt)*r - (+ (/ a tt) - (* (/ b tt) - r))))) - (t - ;; r = 0 so d is very tiny compared to c. - ;; - (/ (+ a (* d (/ b c))) - tt)))) - (robust-subinternal (a b c d) - (declare (vm::double-double-float a b c d)) - (let* ((r (/ d c)) - (tt (+ c (* d r)))) - ;; e is the real part and f is the imaginary part. We - ;; can use internal-compreal for the imaginary part by - ;; noticing that the imaginary part of (a+i*b)/(c+i*d) is - ;; the same as the real part of (b-i*a)/(c+i*d). - (let ((e (internal-compreal a b c d r tt)) - (f (internal-compreal b (- a) c d r tt))) - (values e - f)))) - - (robust-internal (x y) - (declare (type (complex vm::double-double-float) x y)) - (let ((a (realpart x)) - (b (imagpart x)) - (c (realpart y)) - (d (imagpart y))) - (declare (vm::double-double-float a b c d)) - (flet ((maybe-scale (abs-tst a b c d) - (declare (vm::double-double-float a b c d)) - ;; This implements McGehearty's iteration 3 to - ;; handle the case when some values are too big - ;; and should be scaled down. Also if some - ;; values are too tiny, scale them up. - (let ((abs-a (abs a)) - (abs-b (abs b))) - (if (or (> abs-tst +dd-rbig+) - (> abs-a +dd-rbig+) - (> abs-b +dd-rbig+)) - (setf a (* a 0.5d0) - b (* b 0.5d0) - c (* c 0.5d0) - d (* d 0.5d0)) - (if (< abs-tst +dd-rmin2+) - (setf a (* a +dd-rminscal+) - b (* b +dd-rminscal+) - c (* c +dd-rminscal+) - d (* d +dd-rminscal+)) - (if (or (and (< abs-a +dd-rmin+) - (< abs-b +dd-rmax2+) - (< abs-tst +dd-rmax2+)) - (and (< abs-b +dd-rmin+) - (< abs-a +dd-rmax2+) - (< abs-tst +dd-rmax2+))) - (setf a (* a +dd-rminscal+) - b (* b +dd-rminscal+) - c (* c +dd-rminscal+) - d (* d +dd-rminscal+))))) - (values a b c d)))) - (cond - ((<= (abs d) (abs c)) - ;; |d| <= |c|, so we can use robust-subinternal to - ;; perform the division. - (multiple-value-bind (a b c d) - (maybe-scale (abs c) a b c d) - (multiple-value-bind (e f) - (robust-subinternal a b c d) - (complex e f)))) - (t - ;; |d| > |c|. So, instead compute - ;; - ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2) - ;; - ;; Compare this to (a+i*b)/(c+i*d) and we see that - ;; realpart of the former is the same, but the - ;; imagpart of the former is the negative of the - ;; desired division. - (multiple-value-bind (a b c d) - (maybe-scale (abs d) a b c d) - (multiple-value-bind (e f) - (robust-subinternal b a d c) - (complex e (- f)))))))))) - (let* ((a (realpart x)) - (b (imagpart x)) - (c (realpart y)) - (d (imagpart y)) - (ab (max (abs a) (abs b))) - (cd (max (abs c) (abs d))) - (s 1w0)) - (declare (vm::double-double-float s)) - ;; If a or b is big, scale down a and b. - (when (>= ab +dd-rbig+) - (setf x (/ x 2) - s (* s 2))) - ;; If c or d is big, scale down c and d. - (when (>= cd +dd-rbig+) - (setf y (/ y 2) - s (/ s 2))) - ;; If a or b is tiny, scale up a and b. - (when (<= ab (* +dd-rmin+ +dd-2/eps+)) - (setf x (* x +dd-be+) - s (/ s +dd-be+))) - ;; If c or d is tiny, scale up c and d. - (when (<= cd (* +dd-rmin+ +dd-2/eps+)) - (setf y (* y +dd-be+) - s (* s +dd-be+))) - (* s - (robust-internal x y))))) - ;; Smith's algorithm for complex division for (complex single-float). ;; We convert the parts to double-floats before computing the result. (defun cdiv-single-float (x y) @@ -948,7 +797,10 @@ (((complex double-double-float) (foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float))) + (c::cdiv-double-double-float x + (coerce y '(complex double-double-float))) ;; We should do something better for double-double floats. + #+nil (let ((a (realpart x)) (b (imagpart x)) (c (realpart y)) @@ -969,6 +821,9 @@ (((foreach integer ratio single-float double-float double-double-float (complex rational) (complex single-float) (complex double-float)) (complex double-double-float)) + (c::cdiv-double-double-float (coerce x '(complex double-double-float)) + y) + #+nil (let ((a (realpart x)) (b (imagpart x)) (c (realpart y)) ===================================== src/compiler/float-tran-dd.lisp ===================================== @@ -691,4 +691,166 @@ (kernel:double-double-lo a) (kernel:double-double-hi b) (kernel:double-double-lo b))) + +;; An implementation of Baudin and Smith's robust complex division for +;; double-floats. This is a pretty straightforward translation of the +;; original in https://arxiv.org/pdf/1210.4539. +;; +;; This also includes improvements mentioned in +;; https://lpc.events/event/11/contributions/1005/attachments/856/1625/Complex_.... +;; In particular iteration 1 and 3 are added. Iteration 2 and 4 were +;; not added. The test examples from iteration 2 and 4 didn't change +;; with or without changes added. +(let* ((+dd-eps+ (scale-float 1w0 -104)) + (+dd-rmin+ least-positive-normalized-double-double-float) + (+dd-rbig+ (/ most-positive-double-double-float 2)) + (+dd-rmin2+ (scale-float 1w0 -105)) + (+dd-rminscal+ (scale-float 1w0 102)) + (+dd-rmax2+ (* +dd-rbig+ +dd-rmin2+)) + (+dd-be+ (/ 2 (* +dd-eps+ +dd-eps+))) + (+dd-2/eps+ (/ 2 +dd-eps+))) + (declare (double-double-float +dd-eps+ +dd-rmin+ +dd-rbig+ +dd-rmin2+ + +dd-rminscal+ +dd-rmax2+ +dd-be+ +dd-2/eps+)) + (defun cdiv-double-double-float (x y) + (declare (type (complex double-double-float) x y) + (optimize (speed 2) (safety 0))) + (labels + ((internal-compreal (a b c d r tt) + (declare (double-double-float a b c d r tt)) + ;; Compute the real part of the complex division + ;; (a+ib)/(c+id), assuming |c| <= |d|. r = d/c and tt = 1/(c+d*r). + ;; + ;; The realpart is (a*c+b*d)/(c^2+d^2). + ;; + ;; c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r) + ;; + ;; Then + ;; + ;; (a*c+b*d)/(c^2+d^2) = (a*c+b*d)/(c*(c+d*r)) + ;; = (a + b*d/c)/(c+d*r) + ;; = (a + b*r)/(c + d*r). + ;; + ;; Thus tt = (c + d*r). + (cond ((>= (abs r) +dd-rmin+) + (let ((br (* b r))) + (if (/= br 0) + (/ (+ a br) tt) + ;; b*r underflows. Instead, compute + ;; + ;; (a + b*r)*tt = a*tt + b*tt*r + ;; = a*tt + (b*tt)*r + ;; (a + b*r)/tt = a/tt + b/tt*r + ;; = a*tt + (b*tt)*r + (+ (/ a tt) + (* (/ b tt) + r))))) + (t + ;; r = 0 so d is very tiny compared to c. + ;; + (/ (+ a (* d (/ b c))) + tt)))) + (robust-subinternal (a b c d) + (declare (double-double-float a b c d)) + (let* ((r (/ d c)) + (tt (+ c (* d r)))) + ;; e is the real part and f is the imaginary part. We + ;; can use internal-compreal for the imaginary part by + ;; noticing that the imaginary part of (a+i*b)/(c+i*d) is + ;; the same as the real part of (b-i*a)/(c+i*d). + (let ((e (internal-compreal a b c d r tt)) + (f (internal-compreal b (- a) c d r tt))) + (values e + f)))) + + (robust-internal (x y) + (declare (type (complex double-double-float) x y)) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (declare (double-double-float a b c d)) + (flet ((maybe-scale (abs-tst a b c d) + (declare (double-double-float a b c d)) + ;; This implements McGehearty's iteration 3 to + ;; handle the case when some values are too big + ;; and should be scaled down. Also if some + ;; values are too tiny, scale them up. + (let ((abs-a (abs a)) + (abs-b (abs b))) + (if (or (> abs-tst +dd-rbig+) + (> abs-a +dd-rbig+) + (> abs-b +dd-rbig+)) + (setf a (* a 0.5d0) + b (* b 0.5d0) + c (* c 0.5d0) + d (* d 0.5d0)) + (if (< abs-tst +dd-rmin2+) + (setf a (* a +dd-rminscal+) + b (* b +dd-rminscal+) + c (* c +dd-rminscal+) + d (* d +dd-rminscal+)) + (if (or (and (< abs-a +dd-rmin+) + (< abs-b +dd-rmax2+) + (< abs-tst +dd-rmax2+)) + (and (< abs-b +dd-rmin+) + (< abs-a +dd-rmax2+) + (< abs-tst +dd-rmax2+))) + (setf a (* a +dd-rminscal+) + b (* b +dd-rminscal+) + c (* c +dd-rminscal+) + d (* d +dd-rminscal+))))) + (values a b c d)))) + (cond + ((<= (abs d) (abs c)) + ;; |d| <= |c|, so we can use robust-subinternal to + ;; perform the division. + (multiple-value-bind (a b c d) + (maybe-scale (abs c) a b c d) + (multiple-value-bind (e f) + (robust-subinternal a b c d) + (complex e f)))) + (t + ;; |d| > |c|. So, instead compute + ;; + ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2) + ;; + ;; Compare this to (a+i*b)/(c+i*d) and we see that + ;; realpart of the former is the same, but the + ;; imagpart of the former is the negative of the + ;; desired division. + (multiple-value-bind (a b c d) + (maybe-scale (abs d) a b c d) + (multiple-value-bind (e f) + (robust-subinternal b a d c) + (complex e (- f)))))))))) + (let* ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y)) + (ab (max (abs a) (abs b))) + (cd (max (abs c) (abs d))) + (s 1w0)) + (declare (double-double-float s)) + ;; If a or b is big, scale down a and b. + (when (>= ab +dd-rbig+) + (setf x (/ x 2) + s (* s 2))) + ;; If c or d is big, scale down c and d. + (when (>= cd +dd-rbig+) + (setf y (/ y 2) + s (/ s 2))) + ;; If a or b is tiny, scale up a and b. + (when (<= ab (* +dd-rmin+ +dd-2/eps+)) + (setf x (* x +dd-be+) + s (/ s +dd-be+))) + ;; If c or d is tiny, scale up c and d. + (when (<= cd (* +dd-rmin+ +dd-2/eps+)) + (setf y (* y +dd-be+) + s (* s +dd-be+))) + (* s + (robust-internal x y)))))) + +(deftransform / ((x y) ((complex double-double-float) (complex double-double-float)) + *) + `(cdiv-double-double-float x y)) ) ; end progn ===================================== tests/float.lisp ===================================== @@ -509,9 +509,8 @@ (complex (rational (realpart b)) (rational (imagpart b)))) '(complex ext:double-double-float)))) - (let* ((z (kernel::cdiv-double-double-float - (coerce x '(complex ext:double-double-float)) - (coerce y '(complex ext:double-double-float)))) + (let* ((z (/ (coerce x '(complex ext:double-double-float)) + (coerce y '(complex ext:double-double-float)))) (z-true (compute-true x y)) (rel (rel-err z z-true))) (assert-equal expected-rel-w View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/fd4a7a8b0a8d1109e09d898... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/fd4a7a8b0a8d1109e09d898... You're receiving this email because of your account on gitlab.common-lisp.net.