GitLab

at 2026-01-07T20:30:11-08:00

Add cdiv-dd test code

Remove the old test code and add new test code that fits in with the
current scheme used for cdiv-double-float.

     (min (rerr (realpart computed) (realpart expected))

 	 (rerr (imagpart computed) (imagpart expected)))))

-(define-test complex-division.double-double

-  (:tag :issues)

-  (loop for k from 1

-	for test in *test-cases*

-	do

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

-	       test

-	     (declare (ignore expected-rel z-true))

-	     (flet ((compute-true (a b)

-		      ;; Convert a and b to complex rationals, do the

-		      ;; division and convert back to get the true

-		      ;; expected result.

-		      (coerce

-		       (/ (complex (rational (realpart a))

-				   (rational (imagpart a)))

-			  (complex (rational (realpart b))

-				   (rational (imagpart b))))

-		       '(complex ext:double-double-float))))

-	       (let* ((z (/ (coerce x '(complex ext:double-double-float))

-			    (coerce y '(complex ext:double-double-float))))

-		      (z-true (compute-true x y))

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

-		 (assert-equal expected-rel-w

-			       rel

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

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

+

+(defun do-cdiv-dd-test (x y z-true expected-rel)

+  (flet ((compute-true (a b)

+	   ;; Convert a and b to complex rationals, do the

+	   ;; division and convert back to get the true

+	   ;; expected result.

+	   (coerce

+	    (/ (complex (rational (realpart a))

+			(rational (imagpart a)))

+	       (complex (rational (realpart b))

+			(rational (imagpart b))))

+	    '(complex ext:double-double-float))))

+    (let* ((z (/ (coerce x '(complex ext:double-double-float))

+		 (coerce y '(complex ext:double-double-float))))

+	   (z-true (compute-true x y))

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

+      (assert-equal expected-rel

+		    rel

+		    k 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.

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

+    ((frob (name x y z-true rel dd-rel)

        `(define-test ,name

 	  (:tag :issues)

-	  (do-cdiv-test ,x ,y ,z-true ,rel))))

+	  (do-cdiv-test ,x ,y ,z-true ,rel)

+	  (do-cdiv-dd-test ,x ,y ,z-true ,dd-rel))))

   ;; First cases are from Baudin and Smith

   ;; 1

   (frob cdiv.baudin-case.1

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

 	(complex (scale-float 1d0 -1023)

 		 (scale-float -1d0 -1023))

-	53)

+	53

+	106)

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

+	53

+	106)

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

+	53

+	106)

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

+	53

+	106)

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

+	53

+	106)

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

+	53

+	106)

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

+	53

+	106)

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

+	53

+	106)

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

+	53

+	106)

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

+	53

+	106)

   ;; 11

   ;;

   ;; From Maxima.  This was from a (private) email where Maxima used

 	#c(1.2d-311 5.7d-312)

 	#c(3.691993880674614517999740937026568563794896024143749539711267954d301

 	   -1.753697093319947872394996242210428954266103103602859195409591583d301)

-	52)

+	52

+	107)

   ;; 12

   ;;

   ;; Found by ansi tests. z/z should be exactly 1.

 	#c(1.565640716292489d19 0.0d0)

 	#c(1.565640716292489d19 0.0d0)

 	#c(1d0 0)

-	53)

+	53

+	106)

   ;; 13

   ;; Iteration 1.  Without this, we would instead return

   ;;

 		 (parse-hex-float "0x1.98b3fbc1677bbp-697"))

 	(complex (parse-hex-float "0x1.BA8DF8075BCEEp+155")

 		 (parse-hex-float "-0x1.A4AD628DA5B74p-895"))

-	53)

+	53

+	106)

   ;; 14

   ;; Iteration 2.

   (frob cdiv.mcgehearty-iteration.2

 		 (parse-hex-float "0x1.417bd33105808p-256"))

 	(complex (parse-hex-float "0x1.cde593daa4ffep-810")

 		 (parse-hex-float "-0x1.179b9a63df6d3p-799"))

-	52)

+	52

+	106)

   ;; 15

   ;; Iteration 3

   (frob cdiv.mcgehearty-iteration.3

 		 (parse-hex-float "0x1.3d80439e9a119p-731"))

 	(complex (parse-hex-float "-0x1.3b35ed806ae5ap-333")

 		 (parse-hex-float "-0x0.05e01bcbfd9f6p-1022"))

-	53)

+	53

+	106)

   ;; 16

   ;; Iteration 4

   (frob cdiv.mcgehearty-iteration.4

 		 (parse-hex-float "0x1.5bd78c9335899p+1021"))

 	(complex (parse-hex-float "-0x1.423c6ce00c73bp-710")

 		 (parse-hex-float "0x1.d9edcf45bcb0ep-708"))

-	52))

+	52

+	106))

 (define-test complex-division.misc

     (:tag :issue)

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

Commits:

1 changed file:

Changes: