GitLab

at 2025-12-12T06:15:32-08:00

Add reference and comments.

Clean up impl.

Remove the code for computing ulp; I don't want to use that anymore.
We'll use the bits of accuracy used by Baudin and Smith.

Add some comments.

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

+;;

+;; This also includes improvements mentioned in

+;; https://lpc.events/event/11/contributions/1005/attachments/856/1625/Complex_divide.pdf.

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

 (let* (;; This is the value of Scilab's %eps variable.

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

        (rmin least-positive-normalized-double-float)

     (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-%a (string)

   (let* ((sign (if (char= (aref string 0) #\-)

 		   -1

 			   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.

+;; 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.  Tests 13-16 are for examples for improvement

+;; iterations 1-4 from McGehearty.

+;;

+;; Each test is a list of values: x, y, z-true (the value of x/y), and

+;; the bits of accuracy.

 (defparameter *test-cases*

   (list

    ;; 1

 	 52 least-positive-double-float)

    ))

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

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

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

-(defconstant +most-negative-normalized-float-exponent+

-  (nth-value 1 (decode-float least-positive-normalized-double-float)))

-

-(defun power-of-two-p (n)

-  (cond ((integerp n)

-	 (and (> n 0)

-	      (= 0 (logand (abs n) (+ (abs n) -1)))))

-	((floatp n)

-	 (and (> n 0.0)

-	      (= 0.5 (decode-float n))))))

-

-(defun ulp (f)

-  (cond

-    ((= f 0.0)

-     least-positive-double-float)

-    (t

-     (multiple-value-bind (significand expon)

-	 (decode-float f)

-       ;; If the exponent is smaller than the smallest

-       ;; normalized exponent, the ULP is the smallest float.

-       ;; Otherwise, the ULP has an exponent that is

-       ;; float-digits smaller, except when the fraction is a

-       ;; power of two, where we have to increase the exponent

-       ;; by 1.

-       (if (> expon +most-negative-normalized-float-exponent+)

-	   (scale-float 0.5d0

-			(- expon

-			   (float-digits f)

-			   (if (power-of-two-p significand)

-			       0

-			       -1)))

-	   least-positive-double-float)))))

-

-(defun complex-ulp (computed expected)

-  (flet ((real-ulp (f)

-	   (ulp (abs f))))

-    (let ((diff (- computed expected)))

-      (min (real-ulp (realpart diff))

-	   (real-ulp (imagpart diff))))))

-

 (define-test complex-division.double

   (:tag :issues)

   (loop for k from 1

 	       test

 	     (let* ((z (/ x y))

 		    (diff (- z z-true))

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

-		    (ulp (complex-ulp z z-true)))

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

 	       (assert-equal expected-rel

 			     rel

-			     k x y z z-true diff rel)

-	       (assert-true (<= ulp max-ulp)

-			    k x y z z-true diff ulp max-ulp)))))

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

 (define-test complex-division.misc

     (:tag :issue)

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

Commits:

2 changed files:

Changes: