| ... |
... |
@@ -645,144 +645,169 @@ |
|
645
|
645
|
cdiv-double-double-float
|
|
646
|
646
|
two-arg-/))
|
|
647
|
647
|
|
|
648
|
|
-(defun cdiv-double-float (x y)
|
|
649
|
|
- "Accurate division of complex double-float numbers x and y."
|
|
650
|
|
- (declare (type (complex double-float) x y)
|
|
651
|
|
- (optimize (speed 3) (safety 0)))
|
|
652
|
|
- (labels
|
|
653
|
|
- ((internal-compreal (a b c d r tt)
|
|
654
|
|
- (declare (double-float a b c d r tt))
|
|
655
|
|
- ;; Compute the real part of the complex division
|
|
656
|
|
- ;; (a+ib)/(c+id), assuming |c| <= |d|. r = d/c and tt = 1/(c+d*r).
|
|
657
|
|
- ;;
|
|
658
|
|
- ;; The realpart is (a*c+b*d)/(c^2+d^2).
|
|
659
|
|
- ;;
|
|
660
|
|
- ;; c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r)
|
|
661
|
|
- ;;
|
|
662
|
|
- ;; Then
|
|
663
|
|
- ;;
|
|
664
|
|
- ;; (a*c+b*d)/(c^2+d^2) = (a*c+b*d)/(c*(c+d*r))
|
|
665
|
|
- ;; = (a + b*d/c)/(c+d*r)
|
|
666
|
|
- ;; = (a + b*r)/(c + d*r).
|
|
667
|
|
- ;;
|
|
668
|
|
- ;; Thus tt = (c + d*r).
|
|
669
|
|
- (cond ((>= (abs r) +cdiv-rmin+)
|
|
670
|
|
- (let ((br (* b r)))
|
|
671
|
|
- (if (/= br 0)
|
|
672
|
|
- (/ (+ a br) tt)
|
|
673
|
|
- ;; b*r underflows. Instead, compute
|
|
674
|
|
- ;;
|
|
675
|
|
- ;; (a + b*r)*tt = a*tt + b*tt*r
|
|
676
|
|
- ;; = a*tt + (b*tt)*r
|
|
677
|
|
- ;; (a + b*r)/tt = a/tt + b/tt*r
|
|
678
|
|
- ;; = a*tt + (b*tt)*r
|
|
679
|
|
- (+ (/ a tt)
|
|
680
|
|
- (* (/ b tt)
|
|
681
|
|
- r)))))
|
|
682
|
|
- (t
|
|
683
|
|
- ;; r = 0 so d is very tiny compared to c.
|
|
684
|
|
- ;;
|
|
685
|
|
- (/ (+ a (* d (/ b c)))
|
|
686
|
|
- tt))))
|
|
687
|
|
- (robust-subinternal (a b c d)
|
|
688
|
|
- (declare (double-float a b c d))
|
|
689
|
|
- (let* ((r (/ d c))
|
|
690
|
|
- (tt (+ c (* d r))))
|
|
691
|
|
- ;; e is the real part and f is the imaginary part. We
|
|
692
|
|
- ;; can use internal-compreal for the imaginary part by
|
|
693
|
|
- ;; noticing that the imaginary part of (a+i*b)/(c+i*d) is
|
|
694
|
|
- ;; the same as the real part of (b-i*a)/(c+i*d).
|
|
695
|
|
- (let ((e (internal-compreal a b c d r tt))
|
|
696
|
|
- (f (internal-compreal b (- a) c d r tt)))
|
|
697
|
|
- (values e
|
|
698
|
|
- f))))
|
|
699
|
|
- (robust-internal (x y)
|
|
700
|
|
- (declare (type (complex double-float) x y))
|
|
701
|
|
- (let ((a (realpart x))
|
|
702
|
|
- (b (imagpart x))
|
|
703
|
|
- (c (realpart y))
|
|
704
|
|
- (d (imagpart y)))
|
|
705
|
|
- (declare (double-float a b c d))
|
|
706
|
|
- (flet ((maybe-scale (abs-tst a b c d)
|
|
707
|
|
- (declare (double-float a b c d))
|
|
708
|
|
- ;; This implements McGehearty's iteration 3 to
|
|
709
|
|
- ;; handle the case when some values are too big
|
|
710
|
|
- ;; and should be scaled down. Also if some
|
|
711
|
|
- ;; values are too tiny, scale them up.
|
|
712
|
|
- (let ((abs-a (abs a))
|
|
713
|
|
- (abs-b (abs b)))
|
|
714
|
|
- (if (or (> abs-tst +cdiv-rbig+)
|
|
715
|
|
- (> abs-a +cdiv-rbig+)
|
|
716
|
|
- (> abs-b +cdiv-rbig+))
|
|
717
|
|
- (setf a (* a 0.5d0)
|
|
718
|
|
- b (* b 0.5d0)
|
|
719
|
|
- c (* c 0.5d0)
|
|
720
|
|
- d (* d 0.5d0))
|
|
721
|
|
- (if (< abs-tst +cdiv-rmin2+)
|
|
722
|
|
- (setf a (* a +cdiv-rminscal+)
|
|
723
|
|
- b (* b +cdiv-rminscal+)
|
|
724
|
|
- c (* c +cdiv-rminscal+)
|
|
725
|
|
- d (* d +cdiv-rminscal+))
|
|
726
|
|
- (if (or (and (< abs-a +cdiv-rmin+)
|
|
727
|
|
- (< abs-b +cdiv-rmax2+)
|
|
728
|
|
- (< abs-tst +cdiv-rmax2+))
|
|
729
|
|
- (and (< abs-b +cdiv-rmin+)
|
|
730
|
|
- (< abs-a +cdiv-rmax2+)
|
|
731
|
|
- (< abs-tst +cdiv-rmax2+)))
|
|
732
|
|
- (setf a (* a +cdiv-rminscal+)
|
|
733
|
|
- b (* b +cdiv-rminscal+)
|
|
734
|
|
- c (* c +cdiv-rminscal+)
|
|
735
|
|
- d (* d +cdiv-rminscal+)))))
|
|
736
|
|
- (values a b c d))))
|
|
737
|
|
- (cond
|
|
738
|
|
- ((<= (abs d) (abs c))
|
|
739
|
|
- ;; |d| <= |c|, so we can use robust-subinternal to
|
|
740
|
|
- ;; perform the division.
|
|
741
|
|
- (multiple-value-bind (a b c d)
|
|
742
|
|
- (maybe-scale (abs c) a b c d)
|
|
743
|
|
- (multiple-value-bind (e f)
|
|
744
|
|
- (robust-subinternal a b c d)
|
|
745
|
|
- (complex e f))))
|
|
746
|
|
- (t
|
|
747
|
|
- ;; |d| > |c|. So, instead compute
|
|
748
|
|
- ;;
|
|
749
|
|
- ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2)
|
|
750
|
|
- ;;
|
|
751
|
|
- ;; Compare this to (a+i*b)/(c+i*d) and we see that
|
|
752
|
|
- ;; realpart of the former is the same, but the
|
|
753
|
|
- ;; imagpart of the former is the negative of the
|
|
754
|
|
- ;; desired division.
|
|
755
|
|
- (multiple-value-bind (a b c d)
|
|
756
|
|
- (maybe-scale (abs d) a b c d)
|
|
757
|
|
- (multiple-value-bind (e f)
|
|
758
|
|
- (robust-subinternal b a d c)
|
|
759
|
|
- (complex e (- f))))))))))
|
|
760
|
|
- (let* ((a (realpart x))
|
|
761
|
|
- (b (imagpart x))
|
|
762
|
|
- (c (realpart y))
|
|
763
|
|
- (d (imagpart y))
|
|
764
|
|
- (ab (max (abs a) (abs b)))
|
|
765
|
|
- (cd (max (abs c) (abs d)))
|
|
766
|
|
- (s 1d0))
|
|
767
|
|
- (declare (double-float s))
|
|
768
|
|
- ;; If a or b is big, scale down a and b.
|
|
769
|
|
- (when (>= ab +cdiv-rbig+)
|
|
770
|
|
- (setf x (/ x 2)
|
|
771
|
|
- s (* s 2)))
|
|
772
|
|
- ;; If c or d is big, scale down c and d.
|
|
773
|
|
- (when (>= cd +cdiv-rbig+)
|
|
774
|
|
- (setf y (/ y 2)
|
|
775
|
|
- s (/ s 2)))
|
|
776
|
|
- ;; If a or b is tiny, scale up a and b.
|
|
777
|
|
- (when (<= ab (* +cdiv-rmin+ +cdiv-2/eps+))
|
|
778
|
|
- (setf x (* x +cdiv-be+)
|
|
779
|
|
- s (/ s +cdiv-be+)))
|
|
780
|
|
- ;; If c or d is tiny, scale up c and d.
|
|
781
|
|
- (when (<= cd (* +cdiv-rmin+ +cdiv-2/eps+))
|
|
782
|
|
- (setf y (* y +cdiv-be+)
|
|
783
|
|
- s (* s +cdiv-be+)))
|
|
784
|
|
- (* s
|
|
785
|
|
- (robust-internal x y)))))
|
|
|
648
|
+;; The code for both the double-float and double-double-float
|
|
|
649
|
+;; implementation is basically identical except for the constants.
|
|
|
650
|
+;; use a macro to generate both versions.
|
|
|
651
|
+(macrolet
|
|
|
652
|
+ ((frob (type)
|
|
|
653
|
+ (let* ((dd (if (eq type 'double-float)
|
|
|
654
|
+ ""
|
|
|
655
|
+ "DD-"))
|
|
|
656
|
+ (name (symbolicate "CDIV-" type))
|
|
|
657
|
+ (docstring (let ((*print-case* :downcase))
|
|
|
658
|
+ (format nil "Accurate division of complex ~A numbers x and y."
|
|
|
659
|
+ type)))
|
|
|
660
|
+ ;; Create the correct symbols for the constants we need,
|
|
|
661
|
+ ;; as defined above.
|
|
|
662
|
+ (+rmin+ (symbolicate "+CDIV-" dd "RMIN+"))
|
|
|
663
|
+ (+rbig+ (symbolicate "+CDIV-" dd "RBIG+"))
|
|
|
664
|
+ (+rmin2+ (symbolicate "+CDIV-" dd "RMIN2+"))
|
|
|
665
|
+ (+rminscal+ (symbolicate "+CDIV-" dd "RMINSCAL+"))
|
|
|
666
|
+ (+rmax2+ (symbolicate "+CDIV-" dd "RMAX2+"))
|
|
|
667
|
+ (+eps+ (symbolicate "+CDIV-" dd "EPS+"))
|
|
|
668
|
+ (+be+ (symbolicate "+CDIV-" dd "BE+"))
|
|
|
669
|
+ (+2/eps+ (symbolicate "+CDIV-" dd "2/EPS+")))
|
|
|
670
|
+ `(defun ,name (x y)
|
|
|
671
|
+ ,docstring
|
|
|
672
|
+ (declare (type (complex ,type) x y)
|
|
|
673
|
+ (optimize (speed 3) (safety 0)))
|
|
|
674
|
+ (labels
|
|
|
675
|
+ ((internal-compreal (a b c d r tt)
|
|
|
676
|
+ (declare (,type a b c d r tt))
|
|
|
677
|
+ ;; Compute the real part of the complex division
|
|
|
678
|
+ ;; (a+ib)/(c+id), assuming |c| <= |d|. r = d/c and tt = 1/(c+d*r).
|
|
|
679
|
+ ;;
|
|
|
680
|
+ ;; The realpart is (a*c+b*d)/(c^2+d^2).
|
|
|
681
|
+ ;;
|
|
|
682
|
+ ;; c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r)
|
|
|
683
|
+ ;;
|
|
|
684
|
+ ;; Then
|
|
|
685
|
+ ;;
|
|
|
686
|
+ ;; (a*c+b*d)/(c^2+d^2) = (a*c+b*d)/(c*(c+d*r))
|
|
|
687
|
+ ;; = (a + b*d/c)/(c+d*r)
|
|
|
688
|
+ ;; = (a + b*r)/(c + d*r).
|
|
|
689
|
+ ;;
|
|
|
690
|
+ ;; Thus tt = (c + d*r).
|
|
|
691
|
+ (cond ((>= (abs r) ,+rmin+)
|
|
|
692
|
+ (let ((br (* b r)))
|
|
|
693
|
+ (if (/= br 0)
|
|
|
694
|
+ (/ (+ a br) tt)
|
|
|
695
|
+ ;; b*r underflows. Instead, compute
|
|
|
696
|
+ ;;
|
|
|
697
|
+ ;; (a + b*r)*tt = a*tt + b*tt*r
|
|
|
698
|
+ ;; = a*tt + (b*tt)*r
|
|
|
699
|
+ ;; (a + b*r)/tt = a/tt + b/tt*r
|
|
|
700
|
+ ;; = a*tt + (b*tt)*r
|
|
|
701
|
+ (+ (/ a tt)
|
|
|
702
|
+ (* (/ b tt)
|
|
|
703
|
+ r)))))
|
|
|
704
|
+ (t
|
|
|
705
|
+ ;; r = 0 so d is very tiny compared to c.
|
|
|
706
|
+ ;;
|
|
|
707
|
+ (/ (+ a (* d (/ b c)))
|
|
|
708
|
+ tt))))
|
|
|
709
|
+ (robust-subinternal (a b c d)
|
|
|
710
|
+ (declare (,type a b c d))
|
|
|
711
|
+ (let* ((r (/ d c))
|
|
|
712
|
+ (tt (+ c (* d r))))
|
|
|
713
|
+ ;; e is the real part and f is the imaginary part. We
|
|
|
714
|
+ ;; can use internal-compreal for the imaginary part by
|
|
|
715
|
+ ;; noticing that the imaginary part of (a+i*b)/(c+i*d) is
|
|
|
716
|
+ ;; the same as the real part of (b-i*a)/(c+i*d).
|
|
|
717
|
+ (let ((e (internal-compreal a b c d r tt))
|
|
|
718
|
+ (f (internal-compreal b (- a) c d r tt)))
|
|
|
719
|
+ (values e
|
|
|
720
|
+ f))))
|
|
|
721
|
+ (robust-internal (x y)
|
|
|
722
|
+ (declare (type (complex ,type) x y))
|
|
|
723
|
+ (let ((a (realpart x))
|
|
|
724
|
+ (b (imagpart x))
|
|
|
725
|
+ (c (realpart y))
|
|
|
726
|
+ (d (imagpart y)))
|
|
|
727
|
+ (declare (,type a b c d))
|
|
|
728
|
+ (flet ((maybe-scale (abs-tst a b c d)
|
|
|
729
|
+ (declare (,type a b c d))
|
|
|
730
|
+ ;; This implements McGehearty's iteration 3 to
|
|
|
731
|
+ ;; handle the case when some values are too big
|
|
|
732
|
+ ;; and should be scaled down. Also if some
|
|
|
733
|
+ ;; values are too tiny, scale them up.
|
|
|
734
|
+ (let ((abs-a (abs a))
|
|
|
735
|
+ (abs-b (abs b)))
|
|
|
736
|
+ (if (or (> abs-tst ,+rbig+)
|
|
|
737
|
+ (> abs-a ,+rbig+)
|
|
|
738
|
+ (> abs-b ,+rbig+))
|
|
|
739
|
+ (setf a (* a 0.5d0)
|
|
|
740
|
+ b (* b 0.5d0)
|
|
|
741
|
+ c (* c 0.5d0)
|
|
|
742
|
+ d (* d 0.5d0))
|
|
|
743
|
+ (if (< abs-tst ,+rmin2+)
|
|
|
744
|
+ (setf a (* a ,+rminscal+)
|
|
|
745
|
+ b (* b ,+rminscal+)
|
|
|
746
|
+ c (* c ,+rminscal+)
|
|
|
747
|
+ d (* d ,+rminscal+))
|
|
|
748
|
+ (if (or (and (< abs-a ,+rmin+)
|
|
|
749
|
+ (< abs-b ,+rmax2+)
|
|
|
750
|
+ (< abs-tst ,+rmax2+))
|
|
|
751
|
+ (and (< abs-b ,+rmin+)
|
|
|
752
|
+ (< abs-a ,+rmax2+)
|
|
|
753
|
+ (< abs-tst ,+rmax2+)))
|
|
|
754
|
+ (setf a (* a ,+rminscal+)
|
|
|
755
|
+ b (* b ,+rminscal+)
|
|
|
756
|
+ c (* c ,+rminscal+)
|
|
|
757
|
+ d (* d ,+rminscal+)))))
|
|
|
758
|
+ (values a b c d))))
|
|
|
759
|
+ (cond
|
|
|
760
|
+ ((<= (abs d) (abs c))
|
|
|
761
|
+ ;; |d| <= |c|, so we can use robust-subinternal to
|
|
|
762
|
+ ;; perform the division.
|
|
|
763
|
+ (multiple-value-bind (a b c d)
|
|
|
764
|
+ (maybe-scale (abs c) a b c d)
|
|
|
765
|
+ (multiple-value-bind (e f)
|
|
|
766
|
+ (robust-subinternal a b c d)
|
|
|
767
|
+ (complex e f))))
|
|
|
768
|
+ (t
|
|
|
769
|
+ ;; |d| > |c|. So, instead compute
|
|
|
770
|
+ ;;
|
|
|
771
|
+ ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2)
|
|
|
772
|
+ ;;
|
|
|
773
|
+ ;; Compare this to (a+i*b)/(c+i*d) and we see that
|
|
|
774
|
+ ;; realpart of the former is the same, but the
|
|
|
775
|
+ ;; imagpart of the former is the negative of the
|
|
|
776
|
+ ;; desired division.
|
|
|
777
|
+ (multiple-value-bind (a b c d)
|
|
|
778
|
+ (maybe-scale (abs d) a b c d)
|
|
|
779
|
+ (multiple-value-bind (e f)
|
|
|
780
|
+ (robust-subinternal b a d c)
|
|
|
781
|
+ (complex e (- f))))))))))
|
|
|
782
|
+ (let* ((a (realpart x))
|
|
|
783
|
+ (b (imagpart x))
|
|
|
784
|
+ (c (realpart y))
|
|
|
785
|
+ (d (imagpart y))
|
|
|
786
|
+ (ab (max (abs a) (abs b)))
|
|
|
787
|
+ (cd (max (abs c) (abs d)))
|
|
|
788
|
+ (s (coerce 1 ',type)))
|
|
|
789
|
+ (declare (,type s))
|
|
|
790
|
+ ;; If a or b is big, scale down a and b.
|
|
|
791
|
+ (when (>= ab ,+rbig+)
|
|
|
792
|
+ (setf x (/ x 2)
|
|
|
793
|
+ s (* s 2)))
|
|
|
794
|
+ ;; If c or d is big, scale down c and d.
|
|
|
795
|
+ (when (>= cd ,+rbig+)
|
|
|
796
|
+ (setf y (/ y 2)
|
|
|
797
|
+ s (/ s 2)))
|
|
|
798
|
+ ;; If a or b is tiny, scale up a and b.
|
|
|
799
|
+ (when (<= ab (* ,+rmin+ ,+2/eps+))
|
|
|
800
|
+ (setf x (* x ,+be+)
|
|
|
801
|
+ s (/ s ,+be+)))
|
|
|
802
|
+ ;; If c or d is tiny, scale up c and d.
|
|
|
803
|
+ (when (<= cd (* ,+rmin+ ,+2/eps+))
|
|
|
804
|
+ (setf y (* y ,+be+)
|
|
|
805
|
+ s (* s ,+be+)))
|
|
|
806
|
+ (* s
|
|
|
807
|
+ (robust-internal x y))))))))
|
|
|
808
|
+ (frob double-float)
|
|
|
809
|
+ #+double-double
|
|
|
810
|
+ (frob double-double-float))
|
|
786
|
811
|
|
|
787
|
812
|
;; Smith's algorithm for complex division for (complex single-float).
|
|
788
|
813
|
;; We convert the parts to double-floats before computing the result.
|
| ... |
... |
@@ -806,146 +831,6 @@ |
|
806
|
831
|
(f (float (/ (- b (* a r)) denom) 1f0)))
|
|
807
|
832
|
(complex e f))))))
|
|
808
|
833
|
|
|
809
|
|
-(defun cdiv-double-double-float (x y)
|
|
810
|
|
- "Accurate division of complex double-double-float numbers x and y."
|
|
811
|
|
- (declare (type (complex double-double-float) x y)
|
|
812
|
|
- (optimize (speed 2) (safety 0)))
|
|
813
|
|
- (labels
|
|
814
|
|
- ((internal-compreal (a b c d r tt)
|
|
815
|
|
- (declare (double-double-float a b c d r tt))
|
|
816
|
|
- ;; Compute the real part of the complex division
|
|
817
|
|
- ;; (a+ib)/(c+id), assuming |c| <= |d|. r = d/c and tt = 1/(c+d*r).
|
|
818
|
|
- ;;
|
|
819
|
|
- ;; The realpart is (a*c+b*d)/(c^2+d^2).
|
|
820
|
|
- ;;
|
|
821
|
|
- ;; c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r)
|
|
822
|
|
- ;;
|
|
823
|
|
- ;; Then
|
|
824
|
|
- ;;
|
|
825
|
|
- ;; (a*c+b*d)/(c^2+d^2) = (a*c+b*d)/(c*(c+d*r))
|
|
826
|
|
- ;; = (a + b*d/c)/(c+d*r)
|
|
827
|
|
- ;; = (a + b*r)/(c + d*r).
|
|
828
|
|
- ;;
|
|
829
|
|
- ;; Thus tt = (c + d*r).
|
|
830
|
|
- (cond ((>= (abs r) +cdiv-dd-rmin+)
|
|
831
|
|
- (let ((br (* b r)))
|
|
832
|
|
- (if (/= br 0)
|
|
833
|
|
- (/ (+ a br) tt)
|
|
834
|
|
- ;; b*r underflows. Instead, compute
|
|
835
|
|
- ;;
|
|
836
|
|
- ;; (a + b*r)*tt = a*tt + b*tt*r
|
|
837
|
|
- ;; = a*tt + (b*tt)*r
|
|
838
|
|
- ;; (a + b*r)/tt = a/tt + b/tt*r
|
|
839
|
|
- ;; = a*tt + (b*tt)*r
|
|
840
|
|
- (+ (/ a tt)
|
|
841
|
|
- (* (/ b tt)
|
|
842
|
|
- r)))))
|
|
843
|
|
- (t
|
|
844
|
|
- ;; r = 0 so d is very tiny compared to c.
|
|
845
|
|
- ;;
|
|
846
|
|
- (/ (+ a (* d (/ b c)))
|
|
847
|
|
- tt))))
|
|
848
|
|
- (robust-subinternal (a b c d)
|
|
849
|
|
- (declare (double-double-float a b c d))
|
|
850
|
|
- (let* ((r (/ d c))
|
|
851
|
|
- (tt (+ c (* d r))))
|
|
852
|
|
- ;; e is the real part and f is the imaginary part. We
|
|
853
|
|
- ;; can use internal-compreal for the imaginary part by
|
|
854
|
|
- ;; noticing that the imaginary part of (a+i*b)/(c+i*d) is
|
|
855
|
|
- ;; the same as the real part of (b-i*a)/(c+i*d).
|
|
856
|
|
- (let ((e (internal-compreal a b c d r tt))
|
|
857
|
|
- (f (internal-compreal b (- a) c d r tt)))
|
|
858
|
|
- (values e
|
|
859
|
|
- f))))
|
|
860
|
|
-
|
|
861
|
|
- (robust-internal (x y)
|
|
862
|
|
- (declare (type (complex double-double-float) x y))
|
|
863
|
|
- (let ((a (realpart x))
|
|
864
|
|
- (b (imagpart x))
|
|
865
|
|
- (c (realpart y))
|
|
866
|
|
- (d (imagpart y)))
|
|
867
|
|
- (declare (double-double-float a b c d))
|
|
868
|
|
- (flet ((maybe-scale (abs-tst a b c d)
|
|
869
|
|
- (declare (double-double-float a b c d))
|
|
870
|
|
- ;; This implements McGehearty's iteration 3 to
|
|
871
|
|
- ;; handle the case when some values are too big
|
|
872
|
|
- ;; and should be scaled down. Also if some
|
|
873
|
|
- ;; values are too tiny, scale them up.
|
|
874
|
|
- (let ((abs-a (abs a))
|
|
875
|
|
- (abs-b (abs b)))
|
|
876
|
|
- (if (or (> abs-tst +cdiv-dd-rbig+)
|
|
877
|
|
- (> abs-a +cdiv-dd-rbig+)
|
|
878
|
|
- (> abs-b +cdiv-dd-rbig+))
|
|
879
|
|
- (setf a (* a 0.5d0)
|
|
880
|
|
- b (* b 0.5d0)
|
|
881
|
|
- c (* c 0.5d0)
|
|
882
|
|
- d (* d 0.5d0))
|
|
883
|
|
- (if (< abs-tst +cdiv-dd-rmin2+)
|
|
884
|
|
- (setf a (* a +cdiv-dd-rminscal+)
|
|
885
|
|
- b (* b +cdiv-dd-rminscal+)
|
|
886
|
|
- c (* c +cdiv-dd-rminscal+)
|
|
887
|
|
- d (* d +cdiv-dd-rminscal+))
|
|
888
|
|
- (if (or (and (< abs-a +cdiv-dd-rmin+)
|
|
889
|
|
- (< abs-b +cdiv-dd-rmax2+)
|
|
890
|
|
- (< abs-tst +cdiv-dd-rmax2+))
|
|
891
|
|
- (and (< abs-b +cdiv-dd-rmin+)
|
|
892
|
|
- (< abs-a +cdiv-dd-rmax2+)
|
|
893
|
|
- (< abs-tst +cdiv-dd-rmax2+)))
|
|
894
|
|
- (setf a (* a +cdiv-dd-rminscal+)
|
|
895
|
|
- b (* b +cdiv-dd-rminscal+)
|
|
896
|
|
- c (* c +cdiv-dd-rminscal+)
|
|
897
|
|
- d (* d +cdiv-dd-rminscal+)))))
|
|
898
|
|
- (values a b c d))))
|
|
899
|
|
- (cond
|
|
900
|
|
- ((<= (abs d) (abs c))
|
|
901
|
|
- ;; |d| <= |c|, so we can use robust-subinternal to
|
|
902
|
|
- ;; perform the division.
|
|
903
|
|
- (multiple-value-bind (a b c d)
|
|
904
|
|
- (maybe-scale (abs c) a b c d)
|
|
905
|
|
- (multiple-value-bind (e f)
|
|
906
|
|
- (robust-subinternal a b c d)
|
|
907
|
|
- (complex e f))))
|
|
908
|
|
- (t
|
|
909
|
|
- ;; |d| > |c|. So, instead compute
|
|
910
|
|
- ;;
|
|
911
|
|
- ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2)
|
|
912
|
|
- ;;
|
|
913
|
|
- ;; Compare this to (a+i*b)/(c+i*d) and we see that
|
|
914
|
|
- ;; realpart of the former is the same, but the
|
|
915
|
|
- ;; imagpart of the former is the negative of the
|
|
916
|
|
- ;; desired division.
|
|
917
|
|
- (multiple-value-bind (a b c d)
|
|
918
|
|
- (maybe-scale (abs d) a b c d)
|
|
919
|
|
- (multiple-value-bind (e f)
|
|
920
|
|
- (robust-subinternal b a d c)
|
|
921
|
|
- (complex e (- f))))))))))
|
|
922
|
|
- (let* ((a (realpart x))
|
|
923
|
|
- (b (imagpart x))
|
|
924
|
|
- (c (realpart y))
|
|
925
|
|
- (d (imagpart y))
|
|
926
|
|
- (ab (max (abs a) (abs b)))
|
|
927
|
|
- (cd (max (abs c) (abs d)))
|
|
928
|
|
- (s 1w0))
|
|
929
|
|
- (declare (double-double-float s))
|
|
930
|
|
- ;; If a or b is big, scale down a and b.
|
|
931
|
|
- (when (>= ab +cdiv-dd-rbig+)
|
|
932
|
|
- (setf x (/ x 2)
|
|
933
|
|
- s (* s 2)))
|
|
934
|
|
- ;; If c or d is big, scale down c and d.
|
|
935
|
|
- (when (>= cd +cdiv-dd-rbig+)
|
|
936
|
|
- (setf y (/ y 2)
|
|
937
|
|
- s (/ s 2)))
|
|
938
|
|
- ;; If a or b is tiny, scale up a and b.
|
|
939
|
|
- (when (<= ab (* +cdiv-dd-rmin+ +cdiv-dd-2/eps+))
|
|
940
|
|
- (setf x (* x +cdiv-dd-be+)
|
|
941
|
|
- s (/ s +cdiv-dd-be+)))
|
|
942
|
|
- ;; If c or d is tiny, scale up c and d.
|
|
943
|
|
- (when (<= cd (* +cdiv-dd-rmin+ +cdiv-dd-2/eps+))
|
|
944
|
|
- (setf y (* y +cdiv-dd-be+)
|
|
945
|
|
- s (* s +cdiv-dd-be+)))
|
|
946
|
|
- (* s
|
|
947
|
|
- (robust-internal x y)))))
|
|
948
|
|
-
|
|
949
|
834
|
;; Generic implementation of Smith's algorithm.
|
|
950
|
835
|
(defun cdiv-generic (x y)
|
|
951
|
836
|
"Complex division of generic numbers x and y. One of x or y should be
|