-  (:export "CDIV-DOUBLE-DOUBLE-FLOAT")

-  #+double-double

-  (:import-from "C" "CDIV-DOUBLE-DOUBLE-FLOAT")

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

-;;

-;; This is a pretty straightforward change of

-;; kernel::cdiv-double-float for double-double-float.  The constants

-;; may need some tweaking.

-#+nil

-(let* ((+dd-eps+ (scale-float 1w0 -104))

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

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

-       (+dd-rmin2+ (scale-float 1w0 -105))

-       (+dd-rminscal+ (scale-float 1w0 102))

-       (+dd-rmax2+ (* +dd-rbig+ +dd-rmin2+))

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

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

-  (declare (double-double-float +dd-eps+ +dd-rmin+ +dd-rbig+ +dd-rmin2+

-				    +dd-rminscal+ +dd-rmax2+ +dd-be+ +dd-2/eps+))

-  (defun cdiv-double-double-float (x y)

-    (declare (type (complex double-double-float) x y)

-	     (optimize (speed 2) (safety 0)))

-    (labels

-	((internal-compreal (a b c d r tt)

-	   (declare (double-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).

-	   ;;

-	   ;; The realpart is (a*c+b*d)/(c^2+d^2).

-	   ;;

-	   ;;   c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r)

-	   ;;

-	   ;; Then

-	   ;;

-	   ;;   (a*c+b*d)/(c^2+d^2) = (a*c+b*d)/(c*(c+d*r))

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

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

-	   ;;

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

-	   (cond ((>= (abs r) +dd-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.

-		  ;;

-		  (/ (+ a (* d (/ b c)))

-		     tt))))

-	 (robust-subinternal (a b c d)

-	   (declare (double-double-float a b c d))

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

-	   (let ((a (realpart x))

-		 (b (imagpart x))

-		 (c (realpart y))

-		 (d (imagpart y)))

-	     (declare (double-double-float a b c d))

-	     (flet ((maybe-scale (abs-tst a b c d)

-		      (declare (double-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 +dd-rbig+)

-				(> abs-a +dd-rbig+)

-				(> abs-b +dd-rbig+))

-			    (setf a (* a 0.5d0)

-				  b (* b 0.5d0)

-				  c (* c 0.5d0)

-				  d (* d 0.5d0))

-			    (if (< abs-tst +dd-rmin2+)

-				(setf a (* a +dd-rminscal+)

-				      b (* b +dd-rminscal+)

-				      c (* c +dd-rminscal+)

-				      d (* d +dd-rminscal+))

-				(if (or (and (< abs-a +dd-rmin+)

-					     (< abs-b +dd-rmax2+)

-					     (< abs-tst +dd-rmax2+))

-					(and (< abs-b +dd-rmin+)

-					     (< abs-a +dd-rmax2+)

-					     (< abs-tst +dd-rmax2+)))

-				    (setf a (* a +dd-rminscal+)

-					  b (* b +dd-rminscal+)

-					  c (* c +dd-rminscal+)

-					  d (* d +dd-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 1w0))

-	(declare (double-double-float s))

-	;; If a or b is big, scale down a and b.

-	(when (>= ab +dd-rbig+)

-	  (setf x (/ x 2)

-		s (* s 2)))

-	;; If c or d is big, scale down c and d.

-	(when (>= cd +dd-rbig+)

-	  (setf y (/ y 2)

-		s (/ s 2)))

-	;; If a or b is tiny, scale up a and b.

-	(when (<= ab (* +dd-rmin+ +dd-2/eps+))

-	  (setf x (* x +dd-be+)

-		s (/ s +dd-be+)))

-	;; If c or d is tiny, scale up c and d.

-	(when (<= cd (* +dd-rmin+ +dd-2/eps+))

-	  (setf y (* y +dd-be+)

-		s (* s +dd-be+)))

-	(* s

-	   (robust-internal x y))))))

-