| ... |
... |
@@ -774,6 +774,25 @@ |
|
774
|
774
|
(f (float (/ (- b (* a r)) denom) 1f0)))
|
|
775
|
775
|
(complex e f))))))
|
|
776
|
776
|
|
|
|
777
|
+;; Generic implementation of Smith's algorithm.
|
|
|
778
|
+(defun cdiv-generic (x y)
|
|
|
779
|
+ (let ((a (realpart x))
|
|
|
780
|
+ (b (imagpart x))
|
|
|
781
|
+ (c (realpart y))
|
|
|
782
|
+ (d (imagpart y)))
|
|
|
783
|
+ (cond ((< (abs c) (abs d))
|
|
|
784
|
+ (let* ((r (/ c d))
|
|
|
785
|
+ (denom (+ (* c r) d))
|
|
|
786
|
+ (e (/ (+ (* a r) b) denom))
|
|
|
787
|
+ (f (/ (- (* b r) a) denom)))
|
|
|
788
|
+ (canonical-complex e f)))
|
|
|
789
|
+ (t
|
|
|
790
|
+ (let* ((r (/ d c))
|
|
|
791
|
+ (denom (+ c (* d r)))
|
|
|
792
|
+ (e (/ (+ a (* b r)) denom))
|
|
|
793
|
+ (f (/ (- b (* a r)) denom)))
|
|
|
794
|
+ (canonical-complex e f))))))
|
|
|
795
|
+
|
|
777
|
796
|
(defun two-arg-/ (x y)
|
|
778
|
797
|
(number-dispatch ((x number) (y number))
|
|
779
|
798
|
(float-contagion / x y (ratio integer))
|
| ... |
... |
@@ -799,45 +818,17 @@ |
|
799
|
818
|
(foreach (complex rational) (complex single-float) (complex double-float)
|
|
800
|
819
|
(complex double-double-float)))
|
|
801
|
820
|
;; We should do something better for double-double floats.
|
|
802
|
|
- (let ((a (realpart x))
|
|
803
|
|
- (b (imagpart x))
|
|
804
|
|
- (c (realpart y))
|
|
805
|
|
- (d (imagpart y)))
|
|
806
|
|
- (cond ((< (abs c) (abs d))
|
|
807
|
|
- (let* ((r (/ c d))
|
|
808
|
|
- (denom (+ (* c r) d))
|
|
809
|
|
- (e (/ (+ (* a r) b) denom))
|
|
810
|
|
- (f (/ (- (* b r) a) denom)))
|
|
811
|
|
- (canonical-complex e f)))
|
|
812
|
|
- (t
|
|
813
|
|
- (let* ((r (/ d c))
|
|
814
|
|
- (denom (+ c (* d r)))
|
|
815
|
|
- (e (/ (+ a (* b r)) denom))
|
|
816
|
|
- (f (/ (- b (* a r)) denom)))
|
|
817
|
|
- (canonical-complex e f))))))
|
|
|
821
|
+ (cdiv-generic x y))
|
|
818
|
822
|
|
|
819
|
823
|
(((foreach integer ratio single-float double-float double-double-float
|
|
820
|
824
|
(complex rational) (complex single-float) (complex double-float))
|
|
821
|
825
|
(complex double-double-float))
|
|
822
|
|
- (let ((a (realpart x))
|
|
823
|
|
- (b (imagpart x))
|
|
824
|
|
- (c (realpart y))
|
|
825
|
|
- (d (imagpart y)))
|
|
826
|
|
- (cond ((< (abs c) (abs d))
|
|
827
|
|
- (let* ((r (/ c d))
|
|
828
|
|
- (denom (+ (* c r) d))
|
|
829
|
|
- (e (/ (+ (* a r) b) denom))
|
|
830
|
|
- (f (/ (- (* b r) a) denom)))
|
|
831
|
|
- (canonical-complex e f)))
|
|
832
|
|
- (t
|
|
833
|
|
- (let* ((r (/ d c))
|
|
834
|
|
- (denom (+ c (* d r)))
|
|
835
|
|
- (e (/ (+ a (* b r)) denom))
|
|
836
|
|
- (f (/ (- b (* a r)) denom)))
|
|
837
|
|
- (canonical-complex e f))))))
|
|
|
826
|
+ (cdiv-generic x y))
|
|
838
|
827
|
|
|
839
|
828
|
(((foreach integer ratio single-float double-float double-double-float)
|
|
840
|
829
|
(complex rational))
|
|
|
830
|
+ ;; Smith's algorithm, but takes advantage of the fact that the
|
|
|
831
|
+ ;; numerator is a real number and not complex.
|
|
841
|
832
|
(let* ((ry (realpart y))
|
|
842
|
833
|
(iy (imagpart y)))
|
|
843
|
834
|
(if (> (abs ry) (abs iy))
|
| ... |
... |
@@ -853,22 +844,7 @@ |
|
853
|
844
|
(complex rational))
|
|
854
|
845
|
;; We probably don't need to do Smith's algorithm for rationals.
|
|
855
|
846
|
;; A naive implementation of coplex division has no issues.
|
|
856
|
|
- (let ((a (realpart x))
|
|
857
|
|
- (b (imagpart x))
|
|
858
|
|
- (c (realpart y))
|
|
859
|
|
- (d (imagpart y)))
|
|
860
|
|
- (cond ((< (abs c) (abs d))
|
|
861
|
|
- (let* ((r (/ c d))
|
|
862
|
|
- (denom (+ (* c r) d))
|
|
863
|
|
- (e (/ (+ (* a r) b) denom))
|
|
864
|
|
- (f (/ (- (* b r) a) denom)))
|
|
865
|
|
- (canonical-complex e f)))
|
|
866
|
|
- (t
|
|
867
|
|
- (let* ((r (/ d c))
|
|
868
|
|
- (denom (+ c (* d r)))
|
|
869
|
|
- (e (/ (+ a (* b r)) denom))
|
|
870
|
|
- (f (/ (- b (* a r)) denom)))
|
|
871
|
|
- (canonical-complex e f))))))
|
|
|
847
|
+ (cdiv-generic x y))
|
|
872
|
848
|
|
|
873
|
849
|
(((foreach (complex rational) (complex single-float) (complex double-float)
|
|
874
|
850
|
(complex double-double-float))
|