GitLab

at 2026-01-19T13:33:37-08:00

Fix #459:  More accurate complex division for double-double-float

Merge branch 'issue-459-more-accurate-dd-complex-div' into 'master'

Fix #459:  More accurate complex division for double-double-float

Closes #459 and #456

See merge request cmucl/cmucl!338
 ;; 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).

+;;

+;; Make these functions accessible outside the block compilation.

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

+			  cdiv-double-double-float

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

+(eval-when (:compile-toplevel)

+(defmacro define-cdiv (type &key rmin rbig rmin2 rminscal eps)

+  (let* ((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))))

+    `(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 + b*r)/(c + d*r).

 	 ;;

 	 ;; Thus tt = (c + d*r).

-	 (cond ((>= (abs r) +cdiv-rmin+)

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

+		;; r is very small 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))

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

+	   ;; e is the real part and f is the imaginary part.  We can

+	   ;; use 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 ((+rbig+ ,rbig)

+			    (+rmin2+ ,rmin2)

+			    (+rminscal+ ,rminscal)

+			    (+rmax2+ (* ,rbig ,rmin2))

+			    (+rmin+ ,rmin)

+			    (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 (x y)

+	   ,docstring

+	   (declare (type (complex ,type) x y)

+		    ,@opt)

+	   (let ((+rmin+ ,rmin)

+		 (+rbig+ ,rbig)

+		 (+2/eps+ (/ 2 ,eps))

+		 (+be+ (/ 2 (* ,eps ,eps))))

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

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

+	       (multiple-value-bind (e f)

+		   (,internal x y)

+		 (complex (* s e)

+			  (* s f)))))))))

+)

+

+(define-cdiv double-float

+  :rmin least-positive-normalized-double-float

+  :rbig (/ most-positive-double-float 2)

+  :rmin2 (scale-float 1d0 -53)

+  :rminscal (scale-float 1d0 51)

+  ;; This is the value of %eps from Scilab

+  :eps (scale-float 1d0 -52))

+

+;; Same constants but for DOUBLE-DOUBLE-FLOATS.  Some of these aren't

+;; well-defined for DOUBLE-DOUBLE-FLOATS (like eps) so we make our

+;; best guess at what they might be.  Since double-doubles have

+;; twice as many 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 since both

+;; LEAST-POSITIVE-NORMALIZED-DOUBLE-DOUBLE-FLOAT and

+;; MOST-POSITIVE-DOUBLE-DOUBLE-FLOAT have a low component of 0,

+;; 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.  This

+;; also makes some of the operations a bit simpler since operations

+;; between a DOUBLE-DOUBLE-FLOAT and a DOUBLE-FLOAT are simpler than

+;; if both were DOUBLE-DOUBLE-FLOATs.

+(define-cdiv double-double-float

+  :rmin least-positive-normalized-double-float

+  :rbig (/ most-positive-double-float 2)

+  :rmin2 (scale-float 1d0 -106)

+  :rminscal (scale-float 1d0 102)

+  :eps (scale-float 1d0 -104))

+

 ;; Smith's algorithm for complex division for (complex single-float).

 ;; We convert the parts to double-floats before computing the result.

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

+     (cdiv-double-double-float x

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

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

+     (cdiv-double-double-float (coerce x '(complex double-double-float))

+				  y))

     (((foreach integer ratio single-float double-float double-double-float)

       (complex rational))

      (truly-the ,(type-specifier (node-derived-type node))

 		(kernel:%make-double-double-float hi lo))))

+(deftransform / ((a b) ((complex vm::double-double-float) (complex vm::double-double-float))

+		 *)

+  `(kernel::cdiv-double-double-float a b))

+  	      

 (declaim (inline sqr-d))

 (defun sqr-d (a)

   "Square"

             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:

 msgid "Accurate division of complex double-float numbers x and y."

 msgstr ""

+#: src/code/numbers.lisp

+msgid "Accurate division of complex double-double-float numbers x and y."

+msgstr ""

+

 #: src/code/numbers.lisp

 msgid "Accurate division of complex single-float numbers x and y."

 msgstr ""

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

+

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

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

+       `(define-test ,name

+	  (:tag :issues)

+	  (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 1d0)

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

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

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

+	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

+	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

+	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

+	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

+	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

+	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

+	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

+	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

+	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

+	106)

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

+	107)

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

+	106)

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

+	106)

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

+	106)

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

+	106)

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

+	106))

+

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

+

 ;; Issue #458

 (define-test dd-mult-overflow

   (:tag :issues)

Raymond Toy pushed to branch master at cmucl / cmucl

Commits:

5 changed files:

Changes: