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

[Git][cmucl/cmucl] Pushed new branch issue-464-lisp-version-uses-branch-name
by Raymond Toy (＠rtoy) 15 Jan '26

15 Jan '26

Raymond Toy pushed new branch issue-464-lisp-version-uses-branch-name at cmucl / cmucl -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/tree/issue-464-lisp-version-us… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][master] 2 commits: Fix #456: More accurate complex division
by Raymond Toy (＠rtoy) 14 Jan '26

14 Jan '26

Raymond Toy pushed to branch master at cmucl / cmucl Commits: 3416f6d5 by Raymond Toy at 2026-01-08T08:13:59-08:00 Fix #456: More accurate complex division - - - - - f8d90cd0 by Raymond Toy at 2026-01-08T08:13:59-08:00 Merge branch 'issue-456-more-accurate-complex-div' into 'master' Fix #456: More accurate complex division Closes #456 See merge request cmucl/cmucl!335 - - - - - 4 changed files: - src/code/numbers.lisp - src/compiler/float-tran.lisp - src/general-info/release-22a.md - tests/float.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -593,26 +593,250 @@ (build-ratio x y) (build-ratio (truncate x gcd) (truncate y gcd)))))) +;; 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. +(defconstant +cdiv-rmin+ least-positive-normalized-double-float) +(defconstant +cdiv-rbig+ (/ most-positive-double-float 2)) +(defconstant +cdiv-rmin2+ (scale-float 1d0 -53)) +(defconstant +cdiv-rminscal+ (scale-float 1d0 51)) +(defconstant +cdiv-rmax2+ (* +cdiv-rbig+ +cdiv-rmin2+)) +;; This is the value of %eps from Scilab +(defconstant +cdiv-eps+ (scale-float 1d0 -52)) +(defconstant +cdiv-be+ (/ 2 (* +cdiv-eps+ +cdiv-eps+))) +(defconstant +cdiv-2/eps+ (/ 2 +cdiv-eps+)) + +;; Make these functions accessible. cdiv-double-float and +;; cdiv-single-float are used by deftransforms. Of course, two-arg-/ +;; is the interface to division. cdiv-generic isn't used anywhere +;; else. +(declaim (ext:start-block cdiv-double-float cdiv-single-float two-arg-/)) + +(defun cdiv-double-float (x y) + "Accurate division of complex double-float numbers x and 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)) + ;; 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) +cdiv-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-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-float) x y)) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (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 +cdiv-rbig+) + (> abs-a +cdiv-rbig+) + (> abs-b +cdiv-rbig+)) + (setf a (* a 0.5d0) + b (* b 0.5d0) + c (* c 0.5d0) + d (* d 0.5d0)) + (if (< abs-tst +cdiv-rmin2+) + (setf a (* a +cdiv-rminscal+) + b (* b +cdiv-rminscal+) + c (* c +cdiv-rminscal+) + d (* d +cdiv-rminscal+)) + (if (or (and (< abs-a +cdiv-rmin+) + (< abs-b +cdiv-rmax2+) + (< abs-tst +cdiv-rmax2+)) + (and (< abs-b +cdiv-rmin+) + (< abs-a +cdiv-rmax2+) + (< abs-tst +cdiv-rmax2+))) + (setf a (* a +cdiv-rminscal+) + b (* b +cdiv-rminscal+) + c (* c +cdiv-rminscal+) + d (* d +cdiv-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 1d0)) + (declare (double-float s)) + ;; If a or b is big, scale down a and b. + (when (>= ab +cdiv-rbig+) + (setf x (/ x 2) + s (* s 2))) + ;; If c or d is big, scale down c and d. + (when (>= cd +cdiv-rbig+) + (setf y (/ y 2) + s (/ s 2))) + ;; If a or b is tiny, scale up a and b. + (when (<= ab (* +cdiv-rmin+ +cdiv-2/eps+)) + (setf x (* x +cdiv-be+) + s (/ s +cdiv-be+))) + ;; If c or d is tiny, scale up c and d. + (when (<= cd (* +cdiv-rmin+ +cdiv-2/eps+)) + (setf y (* y +cdiv-be+) + s (* s +cdiv-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) + "Accurate division of complex single-float numbers x and y." + (declare (type (complex single-float) x y)) + (let ((a (float (realpart x) 1d0)) + (b (float (imagpart x) 1d0)) + (c (float (realpart y) 1d0)) + (d (float (imagpart y) 1d0))) + (cond ((< (abs c) (abs d)) + (let* ((r (/ c d)) + (denom (+ (* c r) d)) + (e (float (/ (+ (* a r) b) denom) 1f0)) + (f (float (/ (- (* b r) a) denom) 1f0))) + (complex e f))) + (t + (let* ((r (/ d c)) + (denom (+ c (* d r))) + (e (float (/ (+ a (* b r)) denom) 1f0)) + (f (float (/ (- b (* a r)) denom) 1f0))) + (complex e f)))))) + +;; Generic implementation of Smith's algorithm. +(defun cdiv-generic (x y) + "Complex division of generic numbers x and y. One of x or y should be + a complex." + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (cond ((< (abs c) (abs d)) + (let* ((r (/ c d)) + (denom (+ (* c r) d)) + (e (/ (+ (* a r) b) denom)) + (f (/ (- (* b r) a) denom))) + (canonical-complex e f))) + (t + (let* ((r (/ d c)) + (denom (+ c (* d r))) + (e (/ (+ a (* b r)) denom)) + (f (/ (- b (* a r)) denom))) + (canonical-complex e f)))))) (defun two-arg-/ (x y) (number-dispatch ((x number) (y number)) (float-contagion / x y (ratio integer)) - - ((complex complex) - (let* ((rx (realpart x)) - (ix (imagpart x)) - (ry (realpart y)) - (iy (imagpart y))) - (if (> (abs ry) (abs iy)) - (let* ((r (/ iy ry)) - (dn (+ ry (* r iy)))) - (canonical-complex (/ (+ rx (* ix r)) dn) - (/ (- ix (* rx r)) dn))) - (let* ((r (/ ry iy)) - (dn (+ iy (* r ry)))) - (canonical-complex (/ (+ (* rx r) ix) dn) - (/ (- (* ix r) rx) dn)))))) - (((foreach integer ratio single-float double-float) complex) + + (((complex single-float) + (foreach (complex rational) (complex single-float))) + (cdiv-single-float x (coerce y '(complex single-float)))) + (((complex double-float) + (foreach (complex rational) (complex single-float) (complex double-float))) + (cdiv-double-float x (coerce y '(complex double-float)))) + + (((foreach integer ratio single-float (complex rational)) + (complex single-float)) + (cdiv-single-float (coerce x '(complex single-float)) + y)) + + (((foreach integer ratio single-float double-float (complex rational) + (complex single-float)) + (complex double-float)) + (cdiv-double-float (coerce x '(complex double-float)) + y)) + (((complex double-double-float) + (foreach (complex rational) (complex single-float) (complex double-float) + (complex double-double-float))) + ;; We should do something better for double-double floats. + (cdiv-generic x y)) + + (((foreach integer ratio single-float double-float double-double-float + (complex rational) (complex single-float) (complex double-float)) + (complex double-double-float)) + (cdiv-generic x y)) + + (((foreach integer ratio single-float double-float double-double-float) + (complex rational)) + ;; Smith's algorithm, but takes advantage of the fact that the + ;; numerator is a real number and not complex. (let* ((ry (realpart y)) (iy (imagpart y))) (if (> (abs ry) (abs iy)) @@ -624,10 +848,23 @@ (dn (* iy (+ 1 (* r r))))) (canonical-complex (/ (* x r) dn) (/ (- x) dn)))))) - ((complex (or rational float)) + (((complex rational) + (complex rational)) + ;; We probably don't need to do Smith's algorithm for rationals. + ;; A naive implementation of coplex division has no issues. + (cdiv-generic x y)) + + (((foreach (complex rational) (complex single-float) (complex double-float) + (complex double-double-float)) + (or rational float)) (canonical-complex (/ (realpart x) y) (/ (imagpart x) y))) - + + ((double-float + (complex single-float)) + (cdiv-double-float (coerce x '(complex double-float)) + (coerce y '(complex double-float)))) + ((ratio ratio) (let* ((nx (numerator x)) (dx (denominator x)) @@ -656,6 +893,7 @@ (build-ratio (maybe-truncate nx gcd) (* (maybe-truncate y gcd) (denominator x))))))) +(declaim (ext:end-block)) (defun %negate (n) (number-dispatch ((n number)) ===================================== src/compiler/float-tran.lisp ===================================== @@ -1804,46 +1804,13 @@ #+double-double (frob double-double-float)) -#+(and sse2 complex-fp-vops) (macrolet - ((frob (type one) - `(deftransform / ((x y) (,type ,type) * - :policy (> speed space)) - ;; Divide a complex by a complex + ((frob (type name) + `(deftransform / ((x y) ((complex ,type) (complex ,type)) *) + (,name x y)))) + (frob double-float kernel::cdiv-double-float) + (frob single-float kernel::cdiv-single-float)) - ;; Here's how we do a complex division - ;; - ;; Compute (xr + i*xi)/(yr + i*yi) - ;; - ;; Assume |yi| < |yr|. Then - ;; - ;; (xr + i*xi) (xr + i*xi) - ;; ----------- = ----------------- - ;; (yr + i*yi) yr*(1 + i*(yi/yr)) - ;; - ;; (xr + i*xi)*(1 - i*(yi/yr)) - ;; = --------------------------- - ;; yr*(1 + (yi/yr)^2) - ;; - ;; (xr + i*xi)*(1 - i*(yi/yr)) - ;; = --------------------------- - ;; yr + (yi/yr)*yi - ;; - ;; This allows us to use a fast complex multiply followed by - ;; a real division. - '(let* ((ry (realpart y)) - (iy (imagpart y))) - (if (> (abs ry) (abs iy)) - (let* ((r (/ iy ry)) - (dn (+ ry (* r iy)))) - (/ (* x (complex ,one r)) - dn)) - (let* ((r (/ ry iy)) - (dn (+ iy (* r ry)))) - (/ (* x (complex r ,(- one))) - dn))))))) - (frob (complex single-float) 1f0) - (frob (complex double-float) 1d0)) ;;;; Complex contagion: ===================================== src/general-info/release-22a.md ===================================== @@ -34,6 +34,9 @@ public domain. * #446: Use C compiler to get errno values to update UNIX defpackage with errno symbols * #453: Use correct flags for analyzer and always save logs. + * #456: Improve accuracy for division of complex double-floats + using Baudin and Smith's robust complex division algorithm + with improvements by Patrick McGehearty. * #458: Spurious overflow in double-double-float multiply * Other changes: * Improvements to the PCL implementation of CLOS: ===================================== tests/float.lisp ===================================== @@ -348,3 +348,211 @@ (:tag :issues) (assert-equal -2w300 (* -2w300 1w0))) + + + +;; 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-hex-float (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))))) + +;; 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)))) + (if (zerop diff) + (float-digits diff) + (floor (- (log (/ diff (abs e)) 2d0))))))) + (min (rerr (realpart computed) (realpart expected)) + (rerr (imagpart computed) (imagpart expected))))) + +(defun do-cdiv-test (x y z-true expected-rel) + (let* ((z (/ x y)) + (rel (rel-err z z-true))) + (assert-equal expected-rel + rel + x y z z-true rel))) +;; Issue #456: improve accuracy of division of complex double-floats. +;; +;; 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 where (/ z z) didn't produce exactly 1. Tests 13-16 are +;; for examples for improvement iterations 1-4 from McGehearty. +(macrolet + ((frob (name x y z-true rel) + `(define-test ,name + (:tag :issues) + (do-cdiv-test ,x ,y ,z-true ,rel)))) + ;; First cases are from Baudin and Smith + ;; 1 + (frob cdiv.baudin-case.1 + (complex 1d0 1d0) + (complex 1d0 (scale-float 1d0 1023)) + (complex (scale-float 1d0 -1023) + (scale-float -1d0 -1023)) + 53) + ;; 2 + (frob cdiv.baudin-case.2 + (complex 1d0 1d0) + (complex (scale-float 1d0 -1023) (scale-float 1d0 -1023)) + (complex (scale-float 1d0 1023) 0) + 53) + ;; 3 + (frob cdiv.baudin-case.3 + (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) + ;; 4 + (frob cdiv.baudin-case.4.overflow + (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) + (complex 1d0 1d0) + (complex (scale-float 1d0 1023) 0) + 53) + ;; 5 + (frob cdiv.baudin-case.5.underflow-ratio + (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) + ;; 6 + (frob cdiv.baudin-case.6.underflow-realpart + (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) + ;; 7 + (frob cdiv.baudin-case.7.overflow-both-parts + (complex (scale-float 1d0 -347) (scale-float 1d0 -54)) + (complex (scale-float 1d0 -1037) (scale-float 1d0 -1058)) + (complex 3.898125604559113300d289 8.174961907852353577d295) + 53) + ;; 8 + (frob cdiv.baudin-case.8 + (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074)) + (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074)) + (complex 0.6d0 0.2d0) + 53) + ;; 9 + (frob cdiv.baudin-case.9 + (complex (scale-float 1d0 1015) (scale-float 1d0 -989)) + (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) + (complex 0.001953125d0 -0.001953125d0) + 53) + ;; 10 + (frob cdiv.baudin-case.10.improve-imagpart-accuracy + (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) + ;; 11 + ;; + ;; From Maxima. This was from a (private) email where Maxima used + ;; CL:/ to compute the ratio but was not very accurate. + (frob cdiv.maxima-case + #c(5.43d-10 1.13d-100) + #c(1.2d-311 5.7d-312) + #c(3.691993880674614517999740937026568563794896024143749539711267954d301 + -1.753697093319947872394996242210428954266103103602859195409591583d301) + 52) + ;; 12 + ;; + ;; Found by ansi tests. z/z should be exactly 1. + (frob cdiv.ansi-test-z/z + #c(1.565640716292489d19 0.0d0) + #c(1.565640716292489d19 0.0d0) + #c(1d0 0) + 53) + ;; 13 + ;; Iteration 1. Without this, we would instead return + ;; + ;; (complex (parse-hex-float "0x1.ba8df8075bceep+155") + ;; (parse-hex-float "-0x1.a4ad6329485f0p-895")) + ;; + ;; whose imaginary part is quite a bit off. + (frob cdiv.mcgehearty-iteration.1 + (complex (parse-hex-float "0x1.73a3dac1d2f1fp+509") + (parse-hex-float "-0x1.c4dba4ba1ee79p-620")) + (complex (parse-hex-float "0x1.adf526c249cf0p+353") + (parse-hex-float "0x1.98b3fbc1677bbp-697")) + (complex (parse-hex-float "0x1.BA8DF8075BCEEp+155") + (parse-hex-float "-0x1.A4AD628DA5B74p-895")) + 53) + ;; 14 + ;; Iteration 2. + (frob cdiv.mcgehearty-iteration.2 + (complex (parse-hex-float "-0x0.000000008e4f8p-1022") + (parse-hex-float "0x0.0000060366ba7p-1022")) + (complex (parse-hex-float "-0x1.605b467369526p-245") + (parse-hex-float "0x1.417bd33105808p-256")) + (complex (parse-hex-float "0x1.cde593daa4ffep-810") + (parse-hex-float "-0x1.179b9a63df6d3p-799")) + 52) + ;; 15 + ;; Iteration 3 + (frob cdiv.mcgehearty-iteration.3 + (complex (parse-hex-float "0x1.cb27eece7c585p-355 ") + (parse-hex-float "0x0.000000223b8a8p-1022")) + (complex (parse-hex-float "-0x1.74e7ed2b9189fp-22") + (parse-hex-float "0x1.3d80439e9a119p-731")) + (complex (parse-hex-float "-0x1.3b35ed806ae5ap-333") + (parse-hex-float "-0x0.05e01bcbfd9f6p-1022")) + 53) + ;; 16 + ;; Iteration 4 + (frob cdiv.mcgehearty-iteration.4 + (complex (parse-hex-float "-0x1.f5c75c69829f0p-530") + (parse-hex-float "-0x1.e73b1fde6b909p+316")) + (complex (parse-hex-float "-0x1.ff96c3957742bp+1023") + (parse-hex-float "0x1.5bd78c9335899p+1021")) + (complex (parse-hex-float "-0x1.423c6ce00c73bp-710") + (parse-hex-float "0x1.d9edcf45bcb0ep-708")) + 52)) + +(define-test complex-division.misc + (:tag :issue) + (let ((num '(1 + 1/2 + 1.0 + 1d0 + #c(1 2) + #c(1.0 2.0) + #c(1d0 2d0) + #c(1w0 2w0)))) + ;; Try all combinations of divisions of different types. This is + ;; primarily to test that we got all the numeric contagion cases + ;; for division in CL:/. + (dolist (x num) + (dolist (y num) + (assert-true (/ x y) + x y))))) + +(define-test complex-division.single + (:tag :issues) + (let* ((x #c(1 2)) + (y (complex (expt 2 127) (expt 2 127))) + (expected (coerce (/ x y) + '(complex single-float)))) + ;; A naive implementation of complex division would cause an + ;; overflow in computing the denominator. + (assert-equal expected + (/ (coerce x '(complex single-float)) + (coerce y '(complex single-float))) + x + y))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/97824b42c485a2316a2c91… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/97824b42c485a2316a2c91… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-461-cdiv-use-4-doubles] Remove an empty LET
by Raymond Toy (＠rtoy) 14 Jan '26

14 Jan '26

Raymond Toy pushed to branch issue-461-cdiv-use-4-doubles at cmucl / cmucl Commits: eb37d9a9 by Raymond Toy at 2026-01-14T12:14:13-08:00 Remove an empty LET Forgot to remove this empty LET with no bindings. - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -729,59 +729,58 @@ (defun ,internal (a b c d) (declare (,type a b c d) ,@opt) - (let () - (flet ((maybe-scale (abs-tst a b c d) - (declare (,type 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+)) - (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+))))) - (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) - (,subinternal a b c d))) - (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) - (,subinternal b a d c) - (values e (- f))))))))) + (flet ((maybe-scale (abs-tst a b c d) + (declare (,type 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+)) + (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+))))) + (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) + (,subinternal a b c d))) + (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) + (,subinternal b a d c) + (values e (- f)))))))) (defun ,name (a b c d) ,docstring (declare (,type a b c d) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/eb37d9a9b059843f874542c… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/eb37d9a9b059843f874542c… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-461-cdiv-use-4-doubles] 6 commits: Multiply by 0.5 instead of dividing by 2
by Raymond Toy (＠rtoy) 14 Jan '26

14 Jan '26

Raymond Toy pushed to branch issue-461-cdiv-use-4-doubles at cmucl / cmucl Commits: 7fc4e1cb by Raymond Toy at 2026-01-10T19:13:14-08:00 Multiply by 0.5 instead of dividing by 2 For double-floats, we have transforms to convert division by 2 to multiplication by 0.5. No change in precision since this is exact. However, we don't have this for double-double-floats. (We probably should, but that's a different issue.) Instead, just multiply by 0.5d0. This means we can use the faster `mul-dd-d` routine instead of `mul-dd` if we multiplied by 0.5w0. - - - - - abb67a20 by Raymond Toy at 2026-01-11T13:22:31-08:00 Rewrite macro defining the cdiv functions Since we're block-compiling the code, make the macrolet a regular defmacro which then defines 4 functions (compreal, subinternal, internal, and cdiv) that implement the accurate division routine. This makes it easier to debug what's going on. For now, make these new functions visible outside the compilation block. Based on this, we see that the double-float version is nice and only conses when returning the result of the division. The double-double-float cons a lot because it calls all the dd operation code. - - - - - e1a0c380 by Raymond Toy at 2026-01-11T14:31:33-08:00 The variable S can be a double-float As explained in the comment, S can be a double-float in all cases because it's always multiplied by a power of two. Hence for a double-double-float the low part is always zero. - - - - - 23c3c67a by Raymond Toy at 2026-01-11T17:11:18-08:00 Remove old macrolet and don't make internal routines available Don't need the old macrolet anymore so remove it. The internal routines aren't used anywhere else so don't expose them outside the block compilation. - - - - - 3b75ab81 by Raymond Toy at 2026-01-14T08:51:03-08:00 Use double-floats for constants for dd cdiv; update comments All of the constants used for the double-double float version are powers of two and can be exactly represented as doubles instead of double-doubles. This allows slightly faster code we can use simpler routines between a double-double and a double. Add comment about this and also fix up some of the comments that had typos or were incorrect. - - - - - 21115b2d by Raymond Toy at 2026-01-14T12:12:44-08:00 Merge branch 'issue-459-more-accurate-dd-complex-div' into issue-461-cdiv-use-4-doubles Fixed a bunch of merge conflicts too. - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -612,25 +612,33 @@ (defconstant +cdiv-be+ (/ 2 (* +cdiv-eps+ +cdiv-eps+))) (defconstant +cdiv-2/eps+ (/ 2 +cdiv-eps+)) -;; Same constants but for double-double-floats. Some of these aren't -;; well-defined for double-double-floats so we make our best guess at +;; Same constants but for DOUBLE-DOUBLE-FLOATS. Some of these aren't +;; well-defined for DOUBLE-DOUBLE-FLOATS so we make our best guess at ;; what they might be. Since double-doubles have about twice as many -;; bits of precision as a double-float, we generally just double the -;; exponent of the corresponding double-float values above. +;; bits of precision as a DOUBLE-FLOAT, we generally just double the +;; exponent (for SCALE-FLOAT) of the corresponding DOUBLE-FLOAT values +;; above. +;; +;; Note also that both LEAST-POSITIVE-NORMALIZED-DOUBLE-DOUBLE-FLOAT +;; and MOST-POSITIVE-DOUBLE-DOUBLE-FLOAT have a low component of 0, so +;; there's no loss of precision if we use the corresponding +;; DOUBLE-FLOAT values. Likewise, all the other constants here are +;; powers of 2, so the DOUBLE-DOUBLE-FLOAT values can be represented +;; exactly as a DOUBLE-FLOAT. We can use DOUBLE-FLOAT values. (defconstant +cdiv-dd-rmin+ - least-positive-normalized-double-double-float) + least-positive-normalized-double-float) (defconstant +cdiv-dd-rbig+ - (/ most-positive-double-double-float 2)) + (/ most-positive-double-float 2)) (defconstant +cdiv-dd-rmin2+ - (scale-float 1w0 -106)) + (scale-float 1d0 -106)) (defconstant +cdiv-dd-rminscal+ - (scale-float 1w0 102)) + (scale-float 1d0 102)) (defconstant +cdiv-dd-rmax2+ (* +cdiv-dd-rbig+ +cdiv-dd-rmin2+)) ;; Epsilon for double-doubles isn't really well-defined because things ;; like (+ 1w0 1w-200) is a valid double-double float. (defconstant +cdiv-dd-eps+ - (scale-float 1w0 -104)) + (scale-float 1d0 -104)) (defconstant +cdiv-dd-be+ (/ 2 (* +cdiv-dd-eps+ +cdiv-dd-eps+))) (defconstant +cdiv-dd-2/eps+ @@ -648,161 +656,174 @@ ;; The code for both the double-float and double-double-float ;; implementation is basically identical except for the constants. ;; use a macro to generate both versions. -(macrolet - ((frob (type) - (let* ((dd (if (eq type 'double-float) - "" - "DD-")) - (name (symbolicate "CDIV-" type)) - (docstring (let ((*print-case* :downcase)) - (format nil "Accurate division of complex ~A numbers x and y." - type))) - ;; Create the correct symbols for the constants we need, - ;; as defined above. - (+rmin+ (symbolicate "+CDIV-" dd "RMIN+")) - (+rbig+ (symbolicate "+CDIV-" dd "RBIG+")) - (+rmin2+ (symbolicate "+CDIV-" dd "RMIN2+")) - (+rminscal+ (symbolicate "+CDIV-" dd "RMINSCAL+")) - (+rmax2+ (symbolicate "+CDIV-" dd "RMAX2+")) - (+eps+ (symbolicate "+CDIV-" dd "EPS+")) - (+be+ (symbolicate "+CDIV-" dd "BE+")) - (+2/eps+ (symbolicate "+CDIV-" dd "2/EPS+"))) - `(defun ,name (xr xi yr yi) - ,docstring - (declare (,type xr xi yr yi) - (optimize (speed 3) (safety 0))) - (labels - ((internal-compreal (a b c d r tt) - (declare (,type 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) ,+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 (,type 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 (a b c d) - (declare (type ,type a b c d)) - (flet ((maybe-scale (abs-tst a b c d) - (declare (,type 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+)) - (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+))))) - (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* ((xabs (max (abs xr) (abs xi))) - (yabs (max (abs yr) (abs yi))) - (s (coerce 1 ',type))) - (declare (,type s)) - ;; If xr or xi is big, scale down xr and xi. - (when (>= xabs ,+rbig+) - (setf xr (/ xr 2) - xi (/ xi 2) - s (* s 2))) - ;; If yr or yi is big, scale down yr and yi. - (when (>= yabs ,+rbig+) - (setf yr (/ yr 2) - yi (/ yi 2) - s (/ s 2))) - ;; If xr or xi is tiny, scale up xr and xi. - (when (<= xabs (* ,+rmin+ ,+2/eps+)) - (setf xr (* xr ,+be+) - xi (* xi ,+be+) - s (/ s ,+be+))) - ;; If yr or yi is tiny, scale up yr and yi. - (when (<= yabs (* ,+rmin+ ,+2/eps+)) - (setf yr (* yr ,+be+) - yi (* yi ,+be+) - s (* s ,+be+))) - (* s - (robust-internal xr xi yr yi)))))))) - (frob double-float) - #+double-double - (frob double-double-float)) +(eval-when (:compile-toplevel) +(defmacro define-cdiv (type) + (let* ((dd (if (eq type 'double-float) + "" + "DD-")) + (opt '((optimize (speed 3) (safety 0)))) + (name (symbolicate "CDIV-" type)) + (compreal (symbolicate "CDIV-" type "-INTERNAL-COMPREAL")) + (subinternal (symbolicate "CDIV-" type "-ROBUST-SUBINTERNAL")) + (internal (symbolicate "CDIV-" type "-ROBUST-INTERNAL")) + (docstring (let ((*print-case* :downcase)) + (format nil "Accurate division of complex ~A numbers x and y." + type))) + ;; Create the correct symbols for the constants we need, + ;; as defined above. + (+rmin+ (symbolicate "+CDIV-" dd "RMIN+")) + (+rbig+ (symbolicate "+CDIV-" dd "RBIG+")) + (+rmin2+ (symbolicate "+CDIV-" dd "RMIN2+")) + (+rminscal+ (symbolicate "+CDIV-" dd "RMINSCAL+")) + (+rmax2+ (symbolicate "+CDIV-" dd "RMAX2+")) + (+be+ (symbolicate "+CDIV-" dd "BE+")) + (+2/eps+ (symbolicate "+CDIV-" dd "2/EPS+"))) + `(progn + (defun ,compreal (a b c d r tt) + (declare (,type a b c d r tt) + ,@opt) + ;; Computes the real part of the complex division + ;; (a+ib)/(c+id), assuming |d| <= |c|. Then let r = d/c and + ;; tt = (c+d*r). + ;; + ;; The realpart is (a*c+b*d)/(c^2+d^2). + ;; + ;; The denominator is + ;; + ;; 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) ,+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))))) + (t + ;; r is very small so d is very tiny compared to c. + ;; + (/ (+ a (* d (/ b c))) + tt)))) + (defun ,subinternal (a b c d) + (declare (,type a b c d) + ,@opt) + (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 (,compreal a b c d r tt)) + (f (,compreal b (- a) c d r tt))) + (values e f)))) + (defun ,internal (a b c d) + (declare (,type a b c d) + ,@opt) + (let () + (flet ((maybe-scale (abs-tst a b c d) + (declare (,type 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+)) + (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+))))) + (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) + (,subinternal a b c d))) + (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) + (,subinternal b a d c) + (values e (- f))))))))) + (defun ,name (a b c d) + ,docstring + (declare (,type a b c d) + ,@opt) + (let* ((ab (max (abs a) (abs b))) + (cd (max (abs c) (abs d))) + ;; S is always multipled by a power of 2. It's either + ;; 2, 0.5, or +be+, which is also a power of two. + ;; (Either 2^105 or 2^209). Hence, we can just use a + ;; double-float instead of a double-double-float + ;; because the double-double always has 0d0 for the lo + ;; part. + (s 1d0)) + (declare (double-float s)) + ;; If a or b is big, scale down a and b. + (when (>= ab ,+rbig+) + (setf a (* a 0.5d0) + b (* b 0.5d0) + s (* s 2))) + ;; If c or d is big, scale down c and d. + (when (>= cd ,+rbig+) + (setf c (* c 0.5d0) + d (* d 0.5d0) + s (* s 0.5d0))) + ;; If a or b is tiny, scale up a and b. + (when (<= ab (* ,+rmin+ ,+2/eps+)) + (setf a (* a ,+be+) + b (* b ,+be+) + s (/ s ,+be+))) + ;; If c or d is tiny, scale up c and d. + (when (<= cd (* ,+rmin+ ,+2/eps+)) + (setf c (* c ,+be+) + d (* d ,+be+) + s (* s ,+be+))) + (multiple-value-bind (e f) + (,internal a b c d) + (complex (* s e) + (* s f)))))))) +) + +(define-cdiv double-float) +(define-cdiv double-double-float) ;; Smith's algorithm for complex division for (complex single-float). ;; We convert the parts to double-floats before computing the result. View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/66df6947458e1677ffa785… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/66df6947458e1677ffa785… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-459-more-accurate-dd-complex-div] 2 commits: Remove old macrolet and don't make internal routines available
by Raymond Toy (＠rtoy) 14 Jan '26

14 Jan '26

Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: 23c3c67a by Raymond Toy at 2026-01-11T17:11:18-08:00 Remove old macrolet and don't make internal routines available Don't need the old macrolet anymore so remove it. The internal routines aren't used anywhere else so don't expose them outside the block compilation. - - - - - 3b75ab81 by Raymond Toy at 2026-01-14T08:51:03-08:00 Use double-floats for constants for dd cdiv; update comments All of the constants used for the double-double float version are powers of two and can be exactly represented as doubles instead of double-doubles. This allows slightly faster code we can use simpler routines between a double-double and a double. Add comment about this and also fix up some of the comments that had typos or were incorrect. - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -612,25 +612,33 @@ (defconstant +cdiv-be+ (/ 2 (* +cdiv-eps+ +cdiv-eps+))) (defconstant +cdiv-2/eps+ (/ 2 +cdiv-eps+)) -;; Same constants but for double-double-floats. Some of these aren't -;; well-defined for double-double-floats so we make our best guess at +;; Same constants but for DOUBLE-DOUBLE-FLOATS. Some of these aren't +;; well-defined for DOUBLE-DOUBLE-FLOATS so we make our best guess at ;; what they might be. Since double-doubles have about twice as many -;; bits of precision as a double-float, we generally just double the -;; exponent of the corresponding double-float values above. +;; bits of precision as a DOUBLE-FLOAT, we generally just double the +;; exponent (for SCALE-FLOAT) of the corresponding DOUBLE-FLOAT values +;; above. +;; +;; Note also that both LEAST-POSITIVE-NORMALIZED-DOUBLE-DOUBLE-FLOAT +;; and MOST-POSITIVE-DOUBLE-DOUBLE-FLOAT have a low component of 0, so +;; there's no loss of precision if we use the corresponding +;; DOUBLE-FLOAT values. Likewise, all the other constants here are +;; powers of 2, so the DOUBLE-DOUBLE-FLOAT values can be represented +;; exactly as a DOUBLE-FLOAT. We can use DOUBLE-FLOAT values. (defconstant +cdiv-dd-rmin+ - least-positive-normalized-double-double-float) + least-positive-normalized-double-float) (defconstant +cdiv-dd-rbig+ - (/ most-positive-double-double-float 2)) + (/ most-positive-double-float 2)) (defconstant +cdiv-dd-rmin2+ - (scale-float 1w0 -106)) + (scale-float 1d0 -106)) (defconstant +cdiv-dd-rminscal+ - (scale-float 1w0 102)) + (scale-float 1d0 102)) (defconstant +cdiv-dd-rmax2+ (* +cdiv-dd-rbig+ +cdiv-dd-rmin2+)) ;; Epsilon for double-doubles isn't really well-defined because things ;; like (+ 1w0 1w-200) is a valid double-double float. (defconstant +cdiv-dd-eps+ - (scale-float 1w0 -104)) + (scale-float 1d0 -104)) (defconstant +cdiv-dd-be+ (/ 2 (* +cdiv-dd-eps+ +cdiv-dd-eps+))) (defconstant +cdiv-dd-2/eps+ @@ -643,179 +651,11 @@ ;; else. (declaim (ext:start-block cdiv-double-float cdiv-single-float cdiv-double-double-float - two-arg-/ - cdiv-double-float-internal-compreal - cdiv-double-float-robust-subinternal - cdiv-double-float-robust-internal - cdiv-double-double-float-internal-compreal - cdiv-double-double-float-robust-subinternal - cdiv-double-double-float-robust-internal)) + two-arg-/)) ;; The code for both the double-float and double-double-float ;; implementation is basically identical except for the constants. ;; use a macro to generate both versions. -#+nil -(macrolet - ((frob (type) - (let* ((dd (if (eq type 'double-float) - "" - "DD-")) - (name (symbolicate "CDIV-" type)) - (docstring (let ((*print-case* :downcase)) - (format nil "Accurate division of complex ~A numbers x and y." - type))) - ;; Create the correct symbols for the constants we need, - ;; as defined above. - (+rmin+ (symbolicate "+CDIV-" dd "RMIN+")) - (+rbig+ (symbolicate "+CDIV-" dd "RBIG+")) - (+rmin2+ (symbolicate "+CDIV-" dd "RMIN2+")) - (+rminscal+ (symbolicate "+CDIV-" dd "RMINSCAL+")) - (+rmax2+ (symbolicate "+CDIV-" dd "RMAX2+")) - (+eps+ (symbolicate "+CDIV-" dd "EPS+")) - (+be+ (symbolicate "+CDIV-" dd "BE+")) - (+2/eps+ (symbolicate "+CDIV-" dd "2/EPS+"))) - `(defun ,name (x y) - ,docstring - (declare (type (complex ,type) x y) - (optimize (speed 3) (safety 0))) - (labels - ((internal-compreal (a b c d r tt) - (declare (,type 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) ,+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 (,type 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 ,type) x y)) - (let ((a (realpart x)) - (b (imagpart x)) - (c (realpart y)) - (d (imagpart y))) - (declare (,type a b c d)) - (flet ((maybe-scale (abs-tst a b c d) - (declare (,type 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+)) - (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+))))) - (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 (coerce 1 ',type))) - (declare (,type s)) - ;; If a or b is big, scale down a and b. - (when (>= ab ,+rbig+) - (setf x (* x 0.5d0) - s (* s 2))) - ;; If c or d is big, scale down c and d. - (when (>= cd ,+rbig+) - (setf y (* y 0.5d0) - s (* s 0.5d0))) - ;; If a or b is tiny, scale up a and b. - (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+))) - (* s - (robust-internal x y)))))))) - (frob double-float) - #+double-double - (frob double-double-float)) - (eval-when (:compile-toplevel) (defmacro define-cdiv (type) (let* ((dd (if (eq type 'double-float) @@ -842,11 +682,14 @@ (defun ,compreal (a b c d r tt) (declare (,type a b c d r tt) ,@opt) - ;; Compute the real part of the complex division - ;; (a+ib)/(c+id), assuming |c| <= |d|. r = d/c and tt = 1/(c+d*r). + ;; Computes the real part of the complex division + ;; (a+ib)/(c+id), assuming |d| <= |c|. Then let r = d/c and + ;; tt = (c+d*r). ;; ;; The realpart is (a*c+b*d)/(c^2+d^2). ;; + ;; The denominator is + ;; ;; c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r) ;; ;; Then @@ -862,15 +705,12 @@ (/ (+ 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. + ;; r is very small so d is very tiny compared to c. ;; (/ (+ a (* d (/ b c))) tt)))) @@ -931,9 +771,7 @@ ;; perform the division. (multiple-value-bind (a b c d) (maybe-scale (abs c) a b c d) - (multiple-value-bind (e f) - (,subinternal a b c d) - (complex e f)))) + (,subinternal a b c d))) (t ;; |d| > |c|. So, instead compute ;; @@ -947,7 +785,7 @@ (maybe-scale (abs d) a b c d) (multiple-value-bind (e f) (,subinternal b a d c) - (complex e (- f))))))))) + (values e (- f))))))))) (defun ,name (x y) ,docstring (declare (type (complex ,type) x y) @@ -982,14 +820,15 @@ (when (<= cd (* ,+rmin+ ,+2/eps+)) (setf y (* y ,+be+) s (* s ,+be+))) - (* s - (,internal x y))))))) + (multiple-value-bind (e f) + (,internal x y) + (complex (* s e) + (* s f)))))))) ) (define-cdiv double-float) (define-cdiv double-double-float) - ;; 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) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/e1a0c3802ba8a4079d384e… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/e1a0c3802ba8a4079d384e… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-459-more-accurate-dd-complex-div] 2 commits: Rewrite macro defining the cdiv functions
by Raymond Toy (＠rtoy) 12 Jan '26

12 Jan '26

Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: abb67a20 by Raymond Toy at 2026-01-11T13:22:31-08:00 Rewrite macro defining the cdiv functions Since we're block-compiling the code, make the macrolet a regular defmacro which then defines 4 functions (compreal, subinternal, internal, and cdiv) that implement the accurate division routine. This makes it easier to debug what's going on. For now, make these new functions visible outside the compilation block. Based on this, we see that the double-float version is nice and only conses when returning the result of the division. The double-double-float cons a lot because it calls all the dd operation code. - - - - - e1a0c380 by Raymond Toy at 2026-01-11T14:31:33-08:00 The variable S can be a double-float As explained in the comment, S can be a double-float in all cases because it's always multiplied by a power of two. Hence for a double-double-float the low part is always zero. - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -643,11 +643,18 @@ ;; else. (declaim (ext:start-block cdiv-double-float cdiv-single-float cdiv-double-double-float - two-arg-/)) + two-arg-/ + cdiv-double-float-internal-compreal + cdiv-double-float-robust-subinternal + cdiv-double-float-robust-internal + cdiv-double-double-float-internal-compreal + cdiv-double-double-float-robust-subinternal + cdiv-double-double-float-robust-internal)) ;; The code for both the double-float and double-double-float ;; implementation is basically identical except for the constants. ;; use a macro to generate both versions. +#+nil (macrolet ((frob (type) (let* ((dd (if (eq type 'double-float) @@ -809,6 +816,180 @@ #+double-double (frob double-double-float)) +(eval-when (:compile-toplevel) +(defmacro define-cdiv (type) + (let* ((dd (if (eq type 'double-float) + "" + "DD-")) + (opt '((optimize (speed 3) (safety 0)))) + (name (symbolicate "CDIV-" type)) + (compreal (symbolicate "CDIV-" type "-INTERNAL-COMPREAL")) + (subinternal (symbolicate "CDIV-" type "-ROBUST-SUBINTERNAL")) + (internal (symbolicate "CDIV-" type "-ROBUST-INTERNAL")) + (docstring (let ((*print-case* :downcase)) + (format nil "Accurate division of complex ~A numbers x and y." + type))) + ;; Create the correct symbols for the constants we need, + ;; as defined above. + (+rmin+ (symbolicate "+CDIV-" dd "RMIN+")) + (+rbig+ (symbolicate "+CDIV-" dd "RBIG+")) + (+rmin2+ (symbolicate "+CDIV-" dd "RMIN2+")) + (+rminscal+ (symbolicate "+CDIV-" dd "RMINSCAL+")) + (+rmax2+ (symbolicate "+CDIV-" dd "RMAX2+")) + (+be+ (symbolicate "+CDIV-" dd "BE+")) + (+2/eps+ (symbolicate "+CDIV-" dd "2/EPS+"))) + `(progn + (defun ,compreal (a b c d r tt) + (declare (,type a b c d r tt) + ,@opt) + ;; 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) ,+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)))) + (defun ,subinternal (a b c d) + (declare (,type a b c d) + ,@opt) + (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 (,compreal a b c d r tt)) + (f (,compreal b (- a) c d r tt))) + (values e f)))) + (defun ,internal (x y) + (declare (type (complex ,type) x y) + ,@opt) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (declare (,type a b c d)) + (flet ((maybe-scale (abs-tst a b c d) + (declare (,type 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+)) + (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+))))) + (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) + (,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) + (,subinternal b a d c) + (complex e (- f))))))))) + (defun ,name (x y) + ,docstring + (declare (type (complex ,type) x y) + ,@opt) + (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 is always multipled by a power of 2. It's either + ;; 2, 0.5, or +be+, which is also a power of two. + ;; (Either 2^105 or 2^209). Hence, we can just use a + ;; double-float instead of a double-double-float + ;; because the double-double always has 0d0 for the lo + ;; part. + (s 1d0)) + (declare (double-float s)) + ;; If a or b is big, scale down a and b. + (when (>= ab ,+rbig+) + (setf x (* x 0.5d0) + s (* s 2))) + ;; If c or d is big, scale down c and d. + (when (>= cd ,+rbig+) + (setf y (* y 0.5d0) + s (* s 0.5d0))) + ;; If a or b is tiny, scale up a and b. + (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+))) + (* s + (,internal x y))))))) +) + +(define-cdiv double-float) +(define-cdiv double-double-float) + + ;; 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) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/7fc4e1cb7963bc79c697f9… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/7fc4e1cb7963bc79c697f9… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl] Pushed new branch issue-463-double-double-float-comparison
by Raymond Toy (＠rtoy) 11 Jan '26

11 Jan '26

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

1 0

[Git][cmucl/cmucl][issue-459-more-accurate-dd-complex-div] Multiply by 0.5 instead of dividing by 2
by Raymond Toy (＠rtoy) 11 Jan '26

11 Jan '26

Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: 7fc4e1cb by Raymond Toy at 2026-01-10T19:13:14-08:00 Multiply by 0.5 instead of dividing by 2 For double-floats, we have transforms to convert division by 2 to multiplication by 0.5. No change in precision since this is exact. However, we don't have this for double-double-floats. (We probably should, but that's a different issue.) Instead, just multiply by 0.5d0. This means we can use the faster `mul-dd-d` routine instead of `mul-dd` if we multiplied by 0.5w0. - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -789,12 +789,12 @@ (declare (,type s)) ;; If a or b is big, scale down a and b. (when (>= ab ,+rbig+) - (setf x (/ x 2) + (setf x (* x 0.5d0) s (* s 2))) ;; If c or d is big, scale down c and d. (when (>= cd ,+rbig+) - (setf y (/ y 2) - s (/ s 2))) + (setf y (* y 0.5d0) + s (* s 0.5d0))) ;; If a or b is tiny, scale up a and b. (when (<= ab (* ,+rmin+ ,+2/eps+)) (setf x (* x ,+be+) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/7fc4e1cb7963bc79c697f96… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/7fc4e1cb7963bc79c697f96… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl] Pushed new branch issue-461-cdiv-use-4-doubles
by Raymond Toy (＠rtoy) 09 Jan '26

09 Jan '26

Raymond Toy pushed new branch issue-461-cdiv-use-4-doubles at cmucl / cmucl -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/tree/issue-461-cdiv-use-4-doub… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-459-more-accurate-dd-complex-div] Update 22a release notes
by Raymond Toy (＠rtoy) 09 Jan '26

09 Jan '26

Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: f2c487fd by Raymond Toy at 2026-01-09T14:17:07-08:00 Update 22a release notes No CI needed [skip-ci] - - - - - 1 changed file: - src/general-info/release-22a.md Changes: ===================================== src/general-info/release-22a.md ===================================== @@ -38,6 +38,9 @@ public domain. using Baudin and Smith's robust complex division algorithm with improvements by Patrick McGehearty. * #458: Spurious overflow in double-double-float multiply + * #459: Improve accuracy for division of complex + double-double-floats. The same algorithm is used here as + for #456. * Other changes: * Improvements to the PCL implementation of CLOS: * Changes to building procedure: View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/f2c487fdaf1b985fe7c3656… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/f2c487fdaf1b985fe7c3656… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0