Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 949dbaa2 by Raymond Toy at 2025-12-10T20:29:59-08:00 Add improvements from Patrick McGehearty Add suggestion for iteration 1 and add tests for sample cases for the other iterations. This change fixes the test case for iteration 1, which would have failed without this change. For whatever reason, the test from iteration 2 and 4 were already working to at least 52 bits of accuracy. I don't know why. Iteration 3 needs some work. - - - - - 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). - - - - - 2 changed files: - src/code/numbers.lisp - tests/float.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 ((/= r 0) + (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 ===================================== tests/float.lisp ===================================== @@ -343,6 +343,21 @@ (assert-true (typep new-mode 'x86::float-modes)) (assert-equal new-mode (setf (x86::x87-floating-point-modes) new-mode)))) +(defun parse-%a (string) + (let* ((sign (if (char= (aref string 0) #\-) + -1 + 1)) + (dot-posn (position #\. string)) + (p-posn (position #\p string)) + (lead (parse-integer string :start (1- dot-posn) :end dot-posn)) + (frac (parse-integer string :start (1+ dot-posn) :end p-posn :radix 16)) + (exp (parse-integer string :start (1+ p-posn)))) + (* sign + (scale-float (float (+ (ash lead 52) + frac) + 1d0) + (- exp 52))))) + ;; Tests for complex division. From Baudin and Smith, and one test ;; from Maxima. Each test is a list of values: x, y, z-true, ;; bits-of-accuracy, and max-ulp. @@ -412,6 +427,34 @@ #c(1.565640716292489d19 0.0d0) #c(1d0 0) 53 least-positive-double-float) + ;; 13 + ;; Iteration 1. Without this, we would instead return + ;; + ;; (complex (parse-%a "0x1.ba8df8075bceep+155") (parse-%a "-0x1.a4ad6329485f0p-895")) + ;; + ;; whose imaginary part is quite a bit off. + (list (complex (parse-%a "0x1.73a3dac1d2f1fp+509") (parse-%a "-0x1.c4dba4ba1ee79p-620")) + (complex (parse-%a "0x1.adf526c249cf0p+353") (parse-%a "0x1.98b3fbc1677bbp-697")) + (complex (parse-%a "0x1.BA8DF8075BCEEp+155") (parse-%a "-0x1.A4AD628DA5B74p-895")) + 53 least-positive-double-float) + ;; 14 + ;; Iteration 2. + (list (complex (parse-%a "-0x0.000000008e4f8p-1022") (parse-%a "0x0.0000060366ba7p-1022")) + (complex (parse-%a "-0x1.605b467369526p-245") (parse-%a "0x1.417bd33105808p-256")) + (complex (parse-%a "0x1.cde593daa4ffep-810") (parse-%a "-0x1.179b9a63df6d3p-799")) + 52 least-positive-double-float) + ;; 15 + ;; Iteration 3 + (list (complex (parse-%a "0x1.cb27eece7c585p-355 ") (parse-%a "0x0.000000223b8a8p-1022")) + (complex (parse-%a "-0x1.74e7ed2b9189fp-22") (parse-%a "0x1.3d80439e9a119p-731")) + (complex (parse-%a "-0x1.3b35ed806ae5ap-333") (parse-%a "-0x0.05e01bcbfd9f6p-1022")) + 53 least-positive-double-float) + ;; 16 + ;; Iteration 4 + (list (complex (parse-%a "-0x1.f5c75c69829f0p-530") (parse-%a "-0x1.e73b1fde6b909p+316")) + (complex (parse-%a "-0x1.ff96c3957742bp+1023") (parse-%a "0x1.5bd78c9335899p+1021")) + (complex (parse-%a "-0x1.423c6ce00c73bp-710") (parse-%a "0x1.d9edcf45bcb0ep-708")) + 52 least-positive-double-float) )) (defun rel-err (computed expected) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/1aa7c11f6cf9da0be58523e... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/1aa7c11f6cf9da0be58523e... You're receiving this email because of your account on gitlab.common-lisp.net.