| ... |
... |
@@ -596,10 +596,15 @@ |
|
596
|
596
|
;; An implementation of Baudin and Smith's robust complex division for
|
|
597
|
597
|
;; double-floats. This is a pretty straightforward translation of the
|
|
598
|
598
|
;; original in https://arxiv.org/pdf/1210.4539.
|
|
599
|
|
-(let ((ov most-positive-double-float)
|
|
600
|
|
- (un least-positive-normalized-double-float)
|
|
601
|
|
- ;; This is the value of Scilab's %eps variable.
|
|
602
|
|
- (eps (scale-float 1d0 -52)))
|
|
|
599
|
+(let* (;; This is the value of Scilab's %eps variable.
|
|
|
600
|
+ (eps (scale-float 1d0 -52))
|
|
|
601
|
+ (rmin least-positive-normalized-double-float)
|
|
|
602
|
+ (rbig (/ most-positive-double-float 2))
|
|
|
603
|
+ (rmin2 (scale-float 1d0 -53))
|
|
|
604
|
+ (rminscal (scale-float 1d0 51))
|
|
|
605
|
+ (rmax2 (* rbig rmin2))
|
|
|
606
|
+ (be (/ 2 (* eps eps)))
|
|
|
607
|
+ (2/eps (/ 2 eps)))
|
|
603
|
608
|
(defun cdiv-double-float (x y)
|
|
604
|
609
|
(declare (type (complex double-float) x y))
|
|
605
|
610
|
(labels
|
| ... |
... |
@@ -623,7 +628,7 @@ |
|
623
|
628
|
(progn
|
|
624
|
629
|
(format t "a,b,c,d = ~A ~A ~A ~A~%" a b c d)
|
|
625
|
630
|
(format t " r, tt = ~A ~A~%" r tt))
|
|
626
|
|
- (cond ((>= (abs r) un)
|
|
|
631
|
+ (cond ((>= (abs r) rmin)
|
|
627
|
632
|
(let ((br (* b r)))
|
|
628
|
633
|
#+nil
|
|
629
|
634
|
(format t "br = ~A~%" br)
|
| ... |
... |
@@ -660,54 +665,84 @@ |
|
660
|
665
|
(f (internal-compreal b (- a) c d r tt)))
|
|
661
|
666
|
(values e
|
|
662
|
667
|
f))))
|
|
|
668
|
+
|
|
663
|
669
|
(robust-internal (x y)
|
|
664
|
670
|
(declare (type (complex double-float) x y))
|
|
665
|
671
|
(let ((a (realpart x))
|
|
666
|
672
|
(b (imagpart x))
|
|
667
|
673
|
(c (realpart y))
|
|
668
|
674
|
(d (imagpart y)))
|
|
669
|
|
- (cond
|
|
670
|
|
- ((<= (abs d) (abs c))
|
|
671
|
|
- ;; |d| <= |c|, so we can use robust-subinternal to
|
|
672
|
|
- ;; perform the division.
|
|
673
|
|
- (multiple-value-bind (e f)
|
|
674
|
|
- (robust-subinternal a b c d)
|
|
675
|
|
- (complex e f)))
|
|
676
|
|
- (t
|
|
677
|
|
- ;; |d| > |c|. So, instead compute
|
|
678
|
|
- ;;
|
|
679
|
|
- ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2)
|
|
680
|
|
- ;;
|
|
681
|
|
- ;; Compare this to (a+i*b)/(c+i*d) and we see that
|
|
682
|
|
- ;; realpart of the former is the same, but the
|
|
683
|
|
- ;; imagpart of the former is the negative of the
|
|
684
|
|
- ;; desired division.
|
|
685
|
|
- (multiple-value-bind (e f)
|
|
686
|
|
- (robust-subinternal b a d c)
|
|
687
|
|
- (complex e (- f))))))))
|
|
|
675
|
+ (flet ((maybe-scale (abs-tst)
|
|
|
676
|
+ ;; This implements McGehearty's iteration 3 to
|
|
|
677
|
+ ;; handle the case when some values are too big
|
|
|
678
|
+ ;; and should be scaled down. Also if some
|
|
|
679
|
+ ;; values are too tiny, scale them up.
|
|
|
680
|
+ (let ((abs-a (abs a))
|
|
|
681
|
+ (abs-b (abs b)))
|
|
|
682
|
+ (if (or (> abs-tst rbig)
|
|
|
683
|
+ (> abs-a rbig)
|
|
|
684
|
+ (> abs-b rbig))
|
|
|
685
|
+ (setf a (* a 0.5d0)
|
|
|
686
|
+ b (* b 0.5d0)
|
|
|
687
|
+ c (* c 0.5d0)
|
|
|
688
|
+ d (* d 0.5d0))
|
|
|
689
|
+ (if (< abs-tst rmin2)
|
|
|
690
|
+ (setf a (* a rminscal)
|
|
|
691
|
+ b (* b rminscal)
|
|
|
692
|
+ c (* c rminscal)
|
|
|
693
|
+ d (* d rminscal))
|
|
|
694
|
+ (if (or (and (< abs-a rmin)
|
|
|
695
|
+ (< abs-b rmax2)
|
|
|
696
|
+ (< abs-tst rmax2))
|
|
|
697
|
+ (and (< abs-b rmin)
|
|
|
698
|
+ (< abs-a rmax2)
|
|
|
699
|
+ (< abs-tst rmax2)))
|
|
|
700
|
+ (setf a (* a rminscal)
|
|
|
701
|
+ b (* b rminscal)
|
|
|
702
|
+ c (* c rminscal)
|
|
|
703
|
+ d (* d rminscal))))))))
|
|
|
704
|
+ (cond
|
|
|
705
|
+ ((<= (abs d) (abs c))
|
|
|
706
|
+ ;; |d| <= |c|, so we can use robust-subinternal to
|
|
|
707
|
+ ;; perform the division.
|
|
|
708
|
+ (maybe-scale (abs c))
|
|
|
709
|
+ (multiple-value-bind (e f)
|
|
|
710
|
+ (robust-subinternal a b c d)
|
|
|
711
|
+ (complex e f)))
|
|
|
712
|
+ (t
|
|
|
713
|
+ ;; |d| > |c|. So, instead compute
|
|
|
714
|
+ ;;
|
|
|
715
|
+ ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2)
|
|
|
716
|
+ ;;
|
|
|
717
|
+ ;; Compare this to (a+i*b)/(c+i*d) and we see that
|
|
|
718
|
+ ;; realpart of the former is the same, but the
|
|
|
719
|
+ ;; imagpart of the former is the negative of the
|
|
|
720
|
+ ;; desired division.
|
|
|
721
|
+ (maybe-scale (abs d))
|
|
|
722
|
+ (multiple-value-bind (e f)
|
|
|
723
|
+ (robust-subinternal b a d c)
|
|
|
724
|
+ (complex e (- f)))))))))
|
|
688
|
725
|
(let* ((a (realpart x))
|
|
689
|
726
|
(b (imagpart x))
|
|
690
|
727
|
(c (realpart y))
|
|
691
|
728
|
(d (imagpart y))
|
|
692
|
729
|
(ab (max (abs a) (abs b)))
|
|
693
|
730
|
(cd (max (abs c) (abs d)))
|
|
694
|
|
- (bb 2d0)
|
|
695
|
|
- (s 1d0)
|
|
696
|
|
- (be (/ bb (* eps eps))))
|
|
|
731
|
+ (s 1d0))
|
|
697
|
732
|
;; If a or b is big, scale down a and b.
|
|
698
|
|
- (when (>= ab (/ ov 2))
|
|
|
733
|
+ (when (>= ab rbig)
|
|
699
|
734
|
(setf x (/ x 2)
|
|
700
|
735
|
s (* s 2)))
|
|
701
|
736
|
;; If c or d is big, scale down c and d.
|
|
702
|
|
- (when (>= cd (/ ov 2))
|
|
|
737
|
+ (when (>= cd rbig)
|
|
703
|
738
|
(setf y (/ y 2)
|
|
704
|
739
|
s (/ s 2)))
|
|
705
|
740
|
;; If a or b is tiny, scale up a and b.
|
|
706
|
|
- (when (<= ab (* un (/ bb eps)))
|
|
|
741
|
+ (when (<= ab (* rmin 2/eps))
|
|
707
|
742
|
(setf x (* x be)
|
|
708
|
743
|
s (/ s be)))
|
|
709
|
744
|
;; If c or d is tiny, scale up c and d.
|
|
710
|
|
- (when (<= cd (* un (/ bb eps)))
|
|
|
745
|
+ (when (<= cd (* rmin 2/eps))
|
|
711
|
746
|
(setf y (* y be)
|
|
712
|
747
|
s (* s be)))
|
|
713
|
748
|
(* s
|