GitLab

at 2026-01-07T16:30:07-08:00

Add comments and docstrings

Try to give meaningful names to each cdiv test

For each of the tests in `*test-cases*`, break out the test data into
separate individual tests with test names that are (somewhat)
indicative of what's being tested.  This is rather hard, and for some
of Baudin's tests, the article doesn't really seem to say.

Remove old code and cleanup comments

`*test-case*` can be removed now.  Update/reorganize comments
appropriately.

 (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

 ;; 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))

 ;; 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))

 			   1d0)

 		    (- exp 52)))))

-;; 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

-   (list (complex 1d0 1d0)

-	 (complex 1d0 (scale-float 1d0 1023))

-	 (complex (scale-float 1d0 -1023)

-		  (scale-float -1d0 -1023))

-	 53)

-   ;; 2

-   (list (complex 1d0 1d0)

-	 (complex (scale-float 1d0 -1023) (scale-float 1d0 -1023))

-	 (complex (scale-float 1d0 1023) 0)

-	 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)

-   ;; 4

-   (list (complex (scale-float 1d0 1023) (scale-float 1d0 1023))

-	 (complex 1d0 1d0)

-	 (complex (scale-float 1d0 1023) 0)

-	 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)

-   ;; 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)

-   ;; 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)

-   ;; 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)

-   ;; 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)

-   ;; 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)

-   ;; 11

-   ;; From Maxima

-   (list #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.

-   (list #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.

-   (list (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.

-   (list (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

-   (list (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

-   (list (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)

-   ))

-

 ;; 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

     (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.

-(define-test complex-division.double

-  (:tag :issues)

-  (loop for k from 1

-	for test in *test-cases*

-	do

-	   (destructuring-bind (x y z-true expected-rel)

-	       test

-	     (let* ((z (/ x y))

-		    (rel (rel-err z z-true)))

-	       (assert-equal expected-rel

-			     rel

-			     k x y z z-true diff rel)))))

+;;

+;; 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)

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl

Commits:

2 changed files:

Changes: