Raymond Toy pushed to branch issue-456-more-accurate-complex-div-improve at cmucl / cmucl Commits: 43d732df by Raymond Toy at 2025-12-11T08:22:38-08:00 Add iteration 3 from Patrick McGehearty Very straightforward implementation. This could be cleaned up, but I'm making a checkpoint for now. The tests pass. (Well, the ULP test doesn't for test 16 because it's too tight.). - - - - - a75064f2 by Raymond Toy at 2025-12-11T17:01:59-08:00 Refactor iteration 4 code The two sets of code for |c|<=|d| or |d|<=|c| are basically identical except they're testing for c or d. Put that all into a new flet function taking c or d as the arg. - - - - - 11b159dd by Raymond Toy at 2025-12-11T17:10:24-08:00 Add rmin value to match McGehearty code - - - - - 5728d771 by Raymond Toy at 2025-12-11T17:13:25-08:00 Get rid of var ov and un in favor of rbig and rmin - - - - - 8292e387 by Raymond Toy at 2025-12-11T17:16:54-08:00 Some manual common subexpression elimination The arg for `maybe-scale` takes the absolute value so we don't have to keep taking the abs value. Also add local vars to hold the abs value of a and b so we don't have to call abs many times. - - - - - 3d75dda2 by Raymond Toy at 2025-12-11T17:29:33-08:00 Refactor some constants Move be to the list of constants. Get rid of bb which is always 2. This got absorbed into the constants be and 2/eps (new). - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -596,10 +596,15 @@ ;; 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. -(let ((ov most-positive-double-float) - (un least-positive-normalized-double-float) - ;; This is the value of Scilab's %eps variable. - (eps (scale-float 1d0 -52))) +(let* (;; This is the value of Scilab's %eps variable. + (eps (scale-float 1d0 -52)) + (rmin least-positive-normalized-double-float) + (rbig (/ most-positive-double-float 2)) + (rmin2 (scale-float 1d0 -53)) + (rminscal (scale-float 1d0 51)) + (rmax2 (* rbig rmin2)) + (be (/ 2 (* eps eps))) + (2/eps (/ 2 eps))) (defun cdiv-double-float (x y) (declare (type (complex double-float) x y)) (labels @@ -623,7 +628,7 @@ (progn (format t "a,b,c,d = ~A ~A ~A ~A~%" a b c d) (format t " r, tt = ~A ~A~%" r tt)) - (cond ((>= (abs r) un) + (cond ((>= (abs r) rmin) (let ((br (* b r))) #+nil (format t "br = ~A~%" br) @@ -660,54 +665,84 @@ (f (internal-compreal b (- a) c d r tt))) (values e f)))) + (robust-internal (x y) (declare (type (complex double-float) x y)) (let ((a (realpart x)) (b (imagpart x)) (c (realpart y)) (d (imagpart y))) - (cond - ((<= (abs d) (abs c)) - ;; |d| <= |c|, so we can use robust-subinternal to - ;; perform the division. - (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 (e f) - (robust-subinternal b a d c) - (complex e (- f)))))))) + (flet ((maybe-scale (abs-tst) + ;; 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 rbig) + (> abs-a rbig) + (> abs-b rbig)) + (setf a (* a 0.5d0) + b (* b 0.5d0) + c (* c 0.5d0) + d (* d 0.5d0)) + (if (< abs-tst rmin2) + (setf a (* a rminscal) + b (* b rminscal) + c (* c rminscal) + d (* d rminscal)) + (if (or (and (< abs-a rmin) + (< abs-b rmax2) + (< abs-tst rmax2)) + (and (< abs-b rmin) + (< abs-a rmax2) + (< abs-tst rmax2))) + (setf a (* a rminscal) + b (* b rminscal) + c (* c rminscal) + d (* d rminscal)))))))) + (cond + ((<= (abs d) (abs c)) + ;; |d| <= |c|, so we can use robust-subinternal to + ;; perform the division. + (maybe-scale (abs c)) + (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. + (maybe-scale (abs 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))) - (bb 2d0) - (s 1d0) - (be (/ bb (* eps eps)))) + (s 1d0)) ;; If a or b is big, scale down a and b. - (when (>= ab (/ ov 2)) + (when (>= ab rbig) (setf x (/ x 2) s (* s 2))) ;; If c or d is big, scale down c and d. - (when (>= cd (/ ov 2)) + (when (>= cd rbig) (setf y (/ y 2) s (/ s 2))) ;; If a or b is tiny, scale up a and b. - (when (<= ab (* un (/ bb eps))) + (when (<= ab (* rmin 2/eps)) (setf x (* x be) s (/ s be))) ;; If c or d is tiny, scale up c and d. - (when (<= cd (* un (/ bb eps))) + (when (<= cd (* rmin 2/eps)) (setf y (* y be) s (* s be))) (* s View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/949dbaa242c2a77687bcea8... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/949dbaa242c2a77687bcea8... You're receiving this email because of your account on gitlab.common-lisp.net.