GitLab

at 2026-01-19T11:07:59-08:00

Let bind the constants as needed in each routine

This also means the macro needs args to define the values.

Address review comments.

Remove unneeded defconstants and update/clarify some comments

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

-

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

-(defconstant +cdiv-dd-rbig+

-  (/ most-positive-double-float 2))

-(defconstant +cdiv-dd-rmin2+

-  (scale-float 1d0 -106))

-(defconstant +cdiv-dd-rminscal+

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

-(defconstant +cdiv-dd-be+

-  (/ 2 (* +cdiv-dd-eps+ +cdiv-dd-eps+)))

-(defconstant +cdiv-dd-2/eps+

-  (/ 2 +cdiv-dd-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.

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

 ;; implementation is basically identical except for the constants.

 ;; use a macro to generate both versions.

 (eval-when (:compile-toplevel)

-(defmacro define-cdiv (type)

-  (let* ((dd (if (eq type 'double-float)

-		 ""

-		 "DD-"))

-	 (opt '((optimize (speed 3) (safety 0))))

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

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

+			      type))))

     `(progn

        (defun ,compreal (a b c d r tt)

 	 (declare (,type a b c d r tt)

 	 ;;                       = (a + b*r)/(c + d*r).

 	 ;;

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

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

+	 (cond ((>= (abs r) ,rmin)

 		(let ((br (* b r)))

 		  (if (/= br 0)

 		      (/ (+ a br) tt)

 		  ,@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).

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

-	   (multiple-value-bind (e f)

-	       (,internal x y)

-	     (complex (* s e)

-		      (* s f))))))))

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

-(define-cdiv double-double-float)

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

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

Commits:

1 changed file:

Changes: