GitLab

at 2025-12-10T20:29:59-08:00

Add improvements from Patrick McGehearty

Add suggestion for iteration 1 and add tests for sample cases for the
other iterations.  This change fixes the test case for iteration 1,
which would have failed without this change.

For whatever reason, the test from iteration 2 and 4 were already
working to at least 52 bits of accuracy.  I don't know why.

Iteration 3 needs some work.

@@ -596,10 +596,15 @@
  ;; An implementation of Baudin and Smith's robust complex division for
  ;; double-floats.  This is a pretty straightforward translation of the
  ;; original in https://arxiv.org/pdf/1210.4539.
 -(let ((ov most-positive-double-float)
 -      (un least-positive-normalized-double-float)
 -      ;; This is the value of Scilab's %eps variable.
 -      (eps (scale-float 1d0 -52)))
 +(let* (;; This is the value of Scilab's %eps variable.
 +       (eps (scale-float 1d0 -52))
 +       (rmin least-positive-normalized-double-float)
 +       (rbig (/ most-positive-double-float 2))
 +       (rmin2 (scale-float 1d0 -53))
 +       (rminscal (scale-float 1d0 51))
 +       (rmax2 (* rbig rmin2))
 +       (be (/ 2 (* eps eps)))
 +       (2/eps (/ 2 eps)))
    (defun cdiv-double-float (x y)
      (declare (type (complex double-float) x y))
      (labels
@@ -623,7 +628,7 @@
  	   (progn
  	     (format t "a,b,c,d = ~A ~A ~A ~A~%" a b c d)
  	     (format t "  r, tt = ~A ~A~%" r tt))
 -	   (cond ((/= r 0)
 +	   (cond ((>= (abs r) rmin)
  		  (let ((br (* b r)))
  		    #+nil
  		    (format t "br = ~A~%" br)
@@ -660,54 +665,84 @@
  		   (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)))
 -	     (cond
 -	       ((<= (abs d) (abs c))
 -		;; |d| <= |c|, so we can use robust-subinternal to
 -		;; perform the division.
 -		(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 (e f)
 -		    (robust-subinternal b a d c)
 -		  (complex e (- f))))))))
 +	     (flet ((maybe-scale (abs-tst)
 +		      ;; 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))))))))
 +	       (cond
 +		 ((<= (abs d) (abs c))
 +		  ;; |d| <= |c|, so we can use robust-subinternal to
 +		  ;; perform the division.
 +		  (maybe-scale (abs c))
 +		  (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.
 +		  (maybe-scale (abs 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)))
 -	     (bb 2d0)
 -	     (s 1d0)
 -	     (be (/ bb (* eps eps))))
 +	     (s 1d0))
  	;; If a or b is big, scale down a and b.
 -	(when (>= ab (/ ov 2))
 +	(when (>= ab rbig)
  	  (setf x (/ x 2)
  		s (* s 2)))
  	;; If c or d is big, scale down c and d.
 -	(when (>= cd (/ ov 2))
 +	(when (>= cd rbig)
  	  (setf y (/ y 2)
  		s (/ s 2)))
  	;; If a or b is tiny, scale up a and b.
 -	(when (<= ab (* un (/ bb eps)))
 +	(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 (* un (/ bb eps)))
 +	(when (<= cd (* rmin 2/eps))
  	  (setf y (* y be)
  		s (* s be)))
  	(* s

@@ -343,6 +343,21 @@
      (assert-true (typep new-mode 'x86::float-modes))
      (assert-equal new-mode (setf (x86::x87-floating-point-modes) new-mode))))
 +(defun parse-%a (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)))))
++
  ;; Tests for complex division.  From Baudin and Smith, and one test
  ;; from Maxima.  Each test is a list of values: x, y, z-true,
  ;; bits-of-accuracy, and max-ulp.
@@ -412,6 +427,34 @@
  	 #c(1.565640716292489d19 0.0d0)
  	 #c(1d0 0)
 least-positive-double-float)
 +   ;; 13
 +   ;; Iteration 1.  Without this, we would instead return
 +   ;;
 +   ;;   (complex (parse-%a "0x1.ba8df8075bceep+155") (parse-%a "-0x1.a4ad6329485f0p-895"))
 +   ;;
 +   ;; whose imaginary part is quite a bit off.
 +   (list (complex (parse-%a "0x1.73a3dac1d2f1fp+509") (parse-%a "-0x1.c4dba4ba1ee79p-620"))
 +	 (complex (parse-%a "0x1.adf526c249cf0p+353") (parse-%a "0x1.98b3fbc1677bbp-697"))
 +	 (complex (parse-%a "0x1.BA8DF8075BCEEp+155") (parse-%a "-0x1.A4AD628DA5B74p-895"))
 +	 53 least-positive-double-float)
 +   ;; 14
 +   ;; Iteration 2.
 +   (list (complex (parse-%a "-0x0.000000008e4f8p-1022") (parse-%a "0x0.0000060366ba7p-1022"))
 +	 (complex (parse-%a "-0x1.605b467369526p-245") (parse-%a "0x1.417bd33105808p-256"))
 +	 (complex (parse-%a "0x1.cde593daa4ffep-810") (parse-%a "-0x1.179b9a63df6d3p-799"))
 +	 52 least-positive-double-float)
 +   ;; 15
 +   ;; Iteration 3
 +   (list (complex (parse-%a "0x1.cb27eece7c585p-355 ") (parse-%a "0x0.000000223b8a8p-1022"))
 +	 (complex (parse-%a "-0x1.74e7ed2b9189fp-22") (parse-%a "0x1.3d80439e9a119p-731"))
 +	 (complex (parse-%a "-0x1.3b35ed806ae5ap-333") (parse-%a "-0x0.05e01bcbfd9f6p-1022"))
 +	 53 least-positive-double-float)
 +   ;; 16
 +   ;; Iteration 4
 +   (list (complex (parse-%a "-0x1.f5c75c69829f0p-530") (parse-%a "-0x1.e73b1fde6b909p+316"))
 +	 (complex (parse-%a "-0x1.ff96c3957742bp+1023") (parse-%a "0x1.5bd78c9335899p+1021"))
 +	 (complex (parse-%a "-0x1.423c6ce00c73bp-710") (parse-%a "0x1.d9edcf45bcb0ep-708"))
 +	 52 least-positive-double-float)
     ))
  (defun rel-err (computed expected)

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

Commits:

2 changed files:

Changes: