GitLab

at 2026-01-15T09:13:30-08:00

Update cmucl.pot for new docstrings

Forgot to update this when accurate complex division was committed.

[skip-ci]

Fix #463:  Handle double-double-float and double-float comparisons

Merge branch 'issue-463-double-double-float-comparison' into 'master'

Fix #463:  Handle double-double-float and double-float comparisons

Closes #463

See merge request cmucl/cmucl!342
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
Merge branch 'master' into issue-461-cdiv-use-4-doubles

Required fixing merge conflicts with master.

Fix bad merge

Forgot to make internal take 4 args.

 ;; 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 (a b c d)

-	 (declare (,type a b c d)

-		  ,@opt)

+	   (declare (,type a b c d)

+		    ,@opt)

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

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

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

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

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

+	 (defun ,name (xr xi yr yi)

+	   ,docstring

+	   (declare (,type xr xi yr yi)

+		    ,@opt)

+	   (let ((+rmin+ ,rmin)

+		 (+rbig+ ,rbig)

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

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

+	     (let* ((max-x (max (abs xr) (abs xi)))

+		    (max-y (max (abs yr) (abs yi)))

+		    ;; S is always multipled by xr power of 2.  It's

+		    ;; multiplied by either 2, 0.5, or +be+, which is

+		    ;; also xr power of two.  (Either 2^105 or 2^209).

+		    ;; Hence, we can just use xr double-float instead

+		    ;; of xr double-double-float because the

+		    ;; double-double always has 0d0 for the lo part.

+		    (s 1d0))

+	       (declare (double-float s))

+	       ;; If xr or xi is big, scale down xr and xi.

+	       (when (>= max-x +rbig+)

+		 (setf xr (* xr 0.5d0)

+		       xi (* xi 0.5d0)

+		       s (* s 2)))

+	       ;; If yr or yi is big, scale down yr and yi.

+	       (when (>= max-y +rbig+)

+		 (setf yr (* yr 0.5d0)

+		       yi (* yi 0.5d0)

+		       s (* s 0.5d0)))

+	       ;; If xr or xi is tiny, scale up xr and xi.

+	       (when (<= max-x (* +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 (<= max-y (* +rmin+ +2/eps+))

+		 (setf yr (* yr +be+)

+		       yi (* yi +be+)

+		       s (* s +be+)))

+	       (multiple-value-bind (e f)

+		   (,internal xr xi yr yi)

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

 	(kernel:double-double-hi b)

 	(kernel:double-double-lo b)))

+(deftransform = ((a b) (vm::double-double-float double-float) *)

+  `(dd= (kernel:double-double-hi a)

+	(kernel:double-double-lo a)

+	b

+	0d0))

+

+(deftransform = ((a b) (double-float vm::double-double-float) *)

+  `(dd= a

+	0d0

+	(kernel:double-double-hi b)

+	(kernel:double-double-lo b)))

+

 (deftransform < ((a b) (vm::double-double-float vm::double-double-float) *)

   `(dd< (kernel:double-double-hi a)

 	(kernel:double-double-hi b)

 	(kernel:double-double-lo b)))

+(deftransform < ((a b) (vm::double-double-float double-float) *)

+  `(dd< (kernel:double-double-hi a)

+	(kernel:double-double-lo a)

+	b

+	0d0))

+

+(deftransform < ((a b) (double-float vm::double-double-float) *)

+  `(dd< a

+	0d0

+	(kernel:double-double-hi b)

+	(kernel:double-double-lo b)))

+

 (deftransform > ((a b) (vm::double-double-float vm::double-double-float) *)

   `(dd> (kernel:double-double-hi a)

 	(kernel:double-double-lo a)

 	(kernel:double-double-hi b)

 	(kernel:double-double-lo b)))

+

+(deftransform > ((a b) (vm::double-double-float double-float) *)

+  `(dd> (kernel:double-double-hi a)

+	(kernel:double-double-lo a)

+	b

+	0d0))

+

+(deftransform > ((a b) (double-float vm::double-double-float) *)

+  `(dd> a

+	0d0

+	(kernel:double-double-hi b)

+	(kernel:double-double-lo b)))

 ) ; end progn

Raymond Toy pushed to branch issue-461-cdiv-use-4-doubles at cmucl / cmucl

Commits:

2 changed files:

Changes: