December 2025 - cmucl-cvs - mailman3.common-lisp.net

[Git][cmucl/cmucl] Pushed new branch issue-458-double-double-mult-overflow
by Raymond Toy (＠rtoy) 13 Dec '25

13 Dec '25

Raymond Toy pushed new branch issue-458-double-double-mult-overflow at cmucl / cmucl -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/tree/issue-458-double-double-m… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl] Pushed new branch issue-457-delete-directory-signals-errors
by Raymond Toy (＠rtoy) 13 Dec '25

13 Dec '25

Raymond Toy pushed new branch issue-457-delete-directory-signals-errors at cmucl / cmucl -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/tree/issue-457-delete-director… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][master] 2 commits: Fix two important typos in unix-get-username and delete-directory
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch master at cmucl / cmucl Commits: fd2e83a1 by Raymond Toy at 2025-12-12T15:38:35-08:00 Fix two important typos in unix-get-username and delete-directory In `unix-get-username`, we were freeing `name` instead of `result`. In `delete-directory`, we misspelled `recursive` as `recusive`. Fix these typos. - - - - - b9494595 by Raymond Toy at 2025-12-12T15:54:05-08:00 Add some tests for unix-get-username and delete-directory Add a simple test that unix-get-username returns something useful. Add a test for delete-directory that tests it deletes directories as specified, include recursively. - - - - - 4 changed files: - src/code/extensions.lisp - src/code/unix.lisp - tests/pathname.lisp - tests/unix.lisp Changes: ===================================== src/code/extensions.lisp ===================================== @@ -676,7 +676,7 @@ subdirectories. Dirname must be a pathname to a directory. Any NAME or TYPE components in Dirname are ignored." (declare (type pathname dirname)) - (when recusive + (when recursive ;; Find all the files or directories in DIRNAME. (dolist (path (directory (merge-pathnames "*.*" dirname))) ;; If the path is a directory, recursively delete the directory. ===================================== src/code/unix.lisp ===================================== @@ -2564,5 +2564,5 @@ (cast result c-call:c-string) nil) status)) - (free-alien name))))) + (free-alien result))))) ===================================== tests/pathname.lisp ===================================== @@ -142,3 +142,16 @@ "/*.*") :truenamep nil :follow-links nil))) (assert-equal dir-tilde dir-home)))) + +(define-test delete-directory + (let ((dir (ensure-directories-exist "tmp/a/b/c/"))) + ;; Verify that the directories were created. + (assert-equal "tmp/a/b/c/" + dir) + ;; Try to delete the directory. It should fail, which we verify + ;; by noting the directory listing is not empty. + (ext::delete-directory (pathname "tmp/")) + (assert-true (directory "tmp/")) + ;; Now recursively delete the directory. + (ext::delete-directory (pathname "tmp/") :recursive t) + (assert-false (directory "tmp/")))) ===================================== tests/unix.lisp ===================================== @@ -88,4 +88,7 @@ (assert-false result) (assert-true (and (integerp errno) (plusp errno))))) - +(define-test unix-get-username + (let ((uid (unix:unix-getuid))) + (assert-true uid) + (assert-true (unix::unix-get-username uid)))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/ae1bdde7deb88dc3204e62… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/ae1bdde7deb88dc3204e62… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] Get rid of debugging prints not needed anymore
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 6f6a8758 by Raymond Toy at 2025-12-12T09:13:22-08:00 Get rid of debugging prints not needed anymore - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -632,14 +632,8 @@ ;; = (a + b*r)/(c + d*r). ;; ;; Thus tt = (c + d*r). - #+nil - (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) +rmin+) (let ((br (* b r))) - #+nil - (format t "br = ~A~%" br) (if (/= br 0) (/ (+ a br) tt) ;; b*r underflows. Instead, compute @@ -654,11 +648,6 @@ (t ;; r = 0 so d is very tiny compared to c. ;; - ;; (a + b*r)/tt = (a + b*(d/c))/tt - #+nil - (progn - (format t "r = 0~%") - (format t "a*tt = ~A~%" (* a tt))) (/ (+ a (* d (/ b c))) tt)))) (robust-subinternal (a b c d) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/6f6a87580e30627db108944… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/6f6a87580e30627db108944… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] Add declarations and set speed 3 and safety 0
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 272ed882 by Raymond Toy at 2025-12-12T09:08:31-08:00 Add declarations and set speed 3 and safety 0 Set speed 3 and safety 0 for `cdiv-double-float`. This causes lots of compiler warnings, so fix those. The means adding some declarations because the compiler can't figure out the type of the numbers. Also change `maybe-scale` to include the components of the complex numbers. The possibly updated values are returned as multiple values. Update callers appropriately. This gets rid of some extraneous boxing notes. - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -602,17 +602,19 @@ ;; 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* (;; 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))) +(let* ((+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+))) + (declare (double-float +eps+ +rmin+ +rbig+ +rmin2+ + +rminscal+ +rmax2+ +be+ +2/eps+)) (defun cdiv-double-float (x y) - (declare (type (complex double-float) x y)) + (declare (type (complex double-float) x y) + (optimize (speed 3) (safety 0))) (labels ((internal-compreal (a b c d r tt) (declare (double-float a b c d r tt)) @@ -634,7 +636,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) rmin) + (cond ((>= (abs r) +rmin+) (let ((br (* b r))) #+nil (format t "br = ~A~%" br) @@ -678,43 +680,47 @@ (b (imagpart x)) (c (realpart y)) (d (imagpart y))) - (flet ((maybe-scale (abs-tst) + (declare (double-float a b c d)) + (flet ((maybe-scale (abs-tst a b c d) + (declare (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 rbig) - (> abs-a rbig) - (> abs-b rbig)) + (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)))))))) + (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+))))) + (values a b c d)))) (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))) + (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 ;; @@ -724,10 +730,11 @@ ;; 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))))))))) + (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)) @@ -735,22 +742,23 @@ (ab (max (abs a) (abs b))) (cd (max (abs c) (abs d))) (s 1d0)) + (declare (double-float s)) ;; If a or b is big, scale down a and b. - (when (>= ab rbig) + (when (>= ab +rbig+) (setf x (/ x 2) s (* s 2))) ;; If c or d is big, scale down c and d. - (when (>= cd rbig) + (when (>= cd +rbig+) (setf y (/ y 2) s (/ s 2))) ;; If a or b is tiny, scale up a and b. - (when (<= ab (* rmin 2/eps)) - (setf x (* x be) - s (/ s be))) + (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 (* rmin 2/eps)) - (setf y (* y be) - s (* s be))) + (when (<= cd (* +rmin+ +2/eps+)) + (setf y (* y +be+) + s (* s +be+))) (* s (robust-internal x y)))))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/272ed88221456dbd46e1aa4… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/272ed88221456dbd46e1aa4… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] Clean up code more
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 5b00b5d8 by Raymond Toy at 2025-12-12T07:04:17-08:00 Clean up code more Remove least-positive-double-float from *test-cases* because we don't use ulp anymore. Remove unused var `diff` for the same reason. - - - - - 1 changed file: - tests/float.lisp Changes: ===================================== tests/float.lisp ===================================== @@ -375,65 +375,65 @@ (complex 1d0 (scale-float 1d0 1023)) (complex (scale-float 1d0 -1023) (scale-float -1d0 -1023)) - 53 least-positive-double-float) + 53) ;; 2 (list (complex 1d0 1d0) (complex (scale-float 1d0 -1023) (scale-float 1d0 -1023)) (complex (scale-float 1d0 1023) 0) - 53 least-positive-double-float) + 53) ;; 3 (list (complex (scale-float 1d0 1023) (scale-float 1d0 -1023)) (complex (scale-float 1d0 677) (scale-float 1d0 -677)) (complex (scale-float 1d0 346) (scale-float -1d0 -1008)) - 53 least-positive-double-float) + 53) ;; 4 (list (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) (complex 1d0 1d0) (complex (scale-float 1d0 1023) 0) - 53 least-positive-double-float) + 53) ;; 5 (list (complex (scale-float 1d0 1020) (scale-float 1d0 -844)) (complex (scale-float 1d0 656) (scale-float 1d0 -780)) (complex (scale-float 1d0 364) (scale-float -1d0 -1072)) - 53 least-positive-double-float) + 53) ;; 6 (list (complex (scale-float 1d0 -71) (scale-float 1d0 1021)) (complex (scale-float 1d0 1001) (scale-float 1d0 -323)) (complex (scale-float 1d0 -1072) (scale-float 1d0 20)) - 53 least-positive-double-float) + 53) ;; 7 (list (complex (scale-float 1d0 -347) (scale-float 1d0 -54)) (complex (scale-float 1d0 -1037) (scale-float 1d0 -1058)) (complex 3.898125604559113300d289 8.174961907852353577d295) - 53 least-positive-double-float) + 53) ;; 8 (list (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074)) (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074)) (complex 0.6d0 0.2d0) - 53 least-positive-double-float) + 53) ;; 9 (list (complex (scale-float 1d0 1015) (scale-float 1d0 -989)) (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) (complex 0.001953125d0 -0.001953125d0) - 53 least-positive-double-float) + 53) ;; 10 (list (complex (scale-float 1d0 -622) (scale-float 1d0 -1071)) (complex (scale-float 1d0 -343) (scale-float 1d0 -798)) (complex 1.02951151789360578d-84 6.97145987515076231d-220) - 53 least-positive-double-float) + 53) ;; 11 ;; From Maxima (list #c(5.43d-10 1.13d-100) #c(1.2d-311 5.7d-312) #c(3.691993880674614517999740937026568563794896024143749539711267954d301 -1.753697093319947872394996242210428954266103103602859195409591583d301) - 52 least-positive-double-float) + 52) ;; 12 ;; Found by ansi tests. z/z should be exactly 1. (list #c(1.565640716292489d19 0.0d0) #c(1.565640716292489d19 0.0d0) #c(1d0 0) - 53 least-positive-double-float) + 53) ;; 13 ;; Iteration 1. Without this, we would instead return ;; @@ -443,25 +443,25 @@ (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) + 53) ;; 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) + 52) ;; 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) + 53) ;; 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) + 52) )) ;; Relative error in terms of bits of accuracy. This is the @@ -482,12 +482,10 @@ (loop for k from 1 for test in *test-cases* do - (destructuring-bind (x y z-true expected-rel max-ulp) + (destructuring-bind (x y z-true expected-rel) test (let* ((z (/ x y)) - (diff (- z z-true)) (rel (rel-err z z-true))) - (assert-equal expected-rel rel k x y z z-true diff rel))))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/5b00b5d83514f7da9cb65aa… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/5b00b5d83514f7da9cb65aa… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] 2 commits: Add reference and comments.
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 8ca7b5cb by Raymond Toy at 2025-12-12T06:15:32-08:00 Add reference and comments. - - - - - 97e0e433 by Raymond Toy at 2025-12-12T06:44:13-08:00 Clean up impl. Remove the code for computing ulp; I don't want to use that anymore. We'll use the bits of accuracy used by Baudin and Smith. Add some comments. - - - - - 2 changed files: - src/code/numbers.lisp - tests/float.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -596,6 +596,12 @@ ;; 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* (;; This is the value of Scilab's %eps variable. (eps (scale-float 1d0 -52)) (rmin least-positive-normalized-double-float) ===================================== tests/float.lisp ===================================== @@ -343,6 +343,9 @@ (assert-true (typep new-mode 'x86::float-modes)) (assert-equal new-mode (setf (x86::x87-floating-point-modes) new-mode)))) +;; Rudimentary code to read C %a formatted numbers that look like +;; "-0x1.c4dba4ba1ee79p-620". We assume STRING is exactly in this +;; format. No error-checking is done. (defun parse-%a (string) (let* ((sign (if (char= (aref string 0) #\-) -1 @@ -358,9 +361,13 @@ 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. +;; Tests for complex division. Tests 1-10 are from Baudin and Smith. +;; Test 11 is an example from Maxima. Test 12 is an example from the +;; ansi-tests. Tests 13-16 are for examples for improvement +;; iterations 1-4 from McGehearty. +;; +;; Each test is a list of values: x, y, z-true (the value of x/y), and +;; the bits of accuracy. (defparameter *test-cases* (list ;; 1 @@ -457,6 +464,10 @@ 52 least-positive-double-float) )) +;; Relative error in terms of bits of accuracy. This is the +;; definition used by Baudin and Smith. A result of 53 means the two +;; numbers have identical bits. For complex numbers, we use the min +;; of the bits of accuracy of the real and imaginary parts. (defun rel-err (computed expected) (flet ((rerr (c e) (let ((diff (abs (- c e)))) @@ -466,46 +477,6 @@ (min (rerr (realpart computed) (realpart expected)) (rerr (imagpart computed) (imagpart expected))))) -(defconstant +most-negative-normalized-float-exponent+ - (nth-value 1 (decode-float least-positive-normalized-double-float))) - -(defun power-of-two-p (n) - (cond ((integerp n) - (and (> n 0) - (= 0 (logand (abs n) (+ (abs n) -1))))) - ((floatp n) - (and (> n 0.0) - (= 0.5 (decode-float n)))))) - -(defun ulp (f) - (cond - ((= f 0.0) - least-positive-double-float) - (t - (multiple-value-bind (significand expon) - (decode-float f) - ;; If the exponent is smaller than the smallest - ;; normalized exponent, the ULP is the smallest float. - ;; Otherwise, the ULP has an exponent that is - ;; float-digits smaller, except when the fraction is a - ;; power of two, where we have to increase the exponent - ;; by 1. - (if (> expon +most-negative-normalized-float-exponent+) - (scale-float 0.5d0 - (- expon - (float-digits f) - (if (power-of-two-p significand) - 0 - -1))) - least-positive-double-float))))) - -(defun complex-ulp (computed expected) - (flet ((real-ulp (f) - (ulp (abs f)))) - (let ((diff (- computed expected))) - (min (real-ulp (realpart diff)) - (real-ulp (imagpart diff)))))) - (define-test complex-division.double (:tag :issues) (loop for k from 1 @@ -515,14 +486,11 @@ test (let* ((z (/ x y)) (diff (- z z-true)) - (rel (rel-err z z-true)) - (ulp (complex-ulp z z-true))) + (rel (rel-err z z-true))) (assert-equal expected-rel rel - k x y z z-true diff rel) - (assert-true (<= ulp max-ulp) - k x y z z-true diff ulp max-ulp))))) + k x y z z-true diff rel))))) (define-test complex-division.misc (:tag :issue) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3d75dda21fa51f25d6f5cf… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3d75dda21fa51f25d6f5cf… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] 7 commits: Add improvements from Patrick McGehearty
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

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/1aa7c11f6cf9da0be58523… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/1aa7c11f6cf9da0be58523… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div-improve] 6 commits: Add iteration 3 from Patrick McGehearty
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

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/949dbaa242c2a77687bcea… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/949dbaa242c2a77687bcea… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl] Pushed new branch issue-456-more-accurate-complex-div-improve
by Raymond Toy (＠rtoy) 11 Dec '25

11 Dec '25

Raymond Toy pushed new branch issue-456-more-accurate-complex-div-improve at cmucl / cmucl -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/tree/issue-456-more-accurate-c… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0