| ... |
... |
@@ -612,6 +612,31 @@ |
|
612
|
612
|
(defconstant +cdiv-be+ (/ 2 (* +cdiv-eps+ +cdiv-eps+)))
|
|
613
|
613
|
(defconstant +cdiv-2/eps+ (/ 2 +cdiv-eps+))
|
|
614
|
614
|
|
|
|
615
|
+;; Same constants but for double-double-floats. Some of these aren't
|
|
|
616
|
+;; well-defined for double-double-floats so we make our best guess at
|
|
|
617
|
+;; what they might be. Since double-doubles have about twice as many
|
|
|
618
|
+;; bits of precision as a double-float, we generally just double the
|
|
|
619
|
+;; exponent of the corresponding double-float values above.
|
|
|
620
|
+(defconstant +cdiv-dd-rmin+
|
|
|
621
|
+ least-positive-normalized-double-double-float)
|
|
|
622
|
+(defconstant +cdiv-dd-rbig+
|
|
|
623
|
+ (/ most-positive-double-double-float 2))
|
|
|
624
|
+(defconstant +cdiv-dd-rmin2+
|
|
|
625
|
+ (scale-float 1w0 -106))
|
|
|
626
|
+(defconstant +cdiv-dd-rminscal+
|
|
|
627
|
+ (scale-float 1w0 102))
|
|
|
628
|
+(defconstant +cdiv-dd-rmax2+
|
|
|
629
|
+ (* +cdiv-dd-rbig+ +cdiv-dd-rmin2+))
|
|
|
630
|
+;; Epsilon for double-doubles isn't really well-defined because things
|
|
|
631
|
+;; like (+ 1w0 1w-200) is a valid double-double float.
|
|
|
632
|
+(defconstant +cdiv-dd-eps+
|
|
|
633
|
+ (scale-float 1w0 -104))
|
|
|
634
|
+(defconstant +cdiv-dd-be+
|
|
|
635
|
+ (/ 2 (* +cdiv-dd-eps+ +cdiv-dd-eps+)))
|
|
|
636
|
+(defconstant +cdiv-dd-2/eps+
|
|
|
637
|
+ (/ 2 +cdiv-dd-eps+))
|
|
|
638
|
+
|
|
|
639
|
+
|
|
615
|
640
|
;; Make these functions accessible. cdiv-double-float and
|
|
616
|
641
|
;; cdiv-single-float are used by deftransforms. Of course, two-arg-/
|
|
617
|
642
|
;; is the interface to division. cdiv-generic isn't used anywhere
|
| ... |
... |
@@ -781,24 +806,8 @@ |
|
781
|
806
|
(f (float (/ (- b (* a r)) denom) 1f0)))
|
|
782
|
807
|
(complex e f))))))
|
|
783
|
808
|
|
|
784
|
|
-(defconstant +cdiv-dd-eps+
|
|
785
|
|
- (scale-float 1w0 -104))
|
|
786
|
|
-(defconstant +cdiv-dd-rmin+
|
|
787
|
|
- least-positive-normalized-double-double-float)
|
|
788
|
|
-(defconstant +cdiv-dd-rbig+
|
|
789
|
|
- (/ most-positive-double-double-float 2))
|
|
790
|
|
-(defconstant +cdiv-dd-rmin2+
|
|
791
|
|
- (scale-float 1w0 -105))
|
|
792
|
|
-(defconstant +cdiv-dd-rminscal+
|
|
793
|
|
- (scale-float 1w0 102))
|
|
794
|
|
-(defconstant +cdiv-dd-rmax2+
|
|
795
|
|
- (* +cdiv-dd-rbig+ +cdiv-dd-rmin2+))
|
|
796
|
|
-(defconstant +cdiv-dd-be+
|
|
797
|
|
- (/ 2 (* +cdiv-dd-eps+ +cdiv-dd-eps+)))
|
|
798
|
|
-(defconstant +cdiv-dd-2/eps+
|
|
799
|
|
- (/ 2 +cdiv-dd-eps+))
|
|
800
|
|
-
|
|
801
|
809
|
(defun cdiv-double-double-float (x y)
|
|
|
810
|
+ "Accurate division of complex double-double-float numbers x and y."
|
|
802
|
811
|
(declare (type (complex double-double-float) x y)
|
|
803
|
812
|
(optimize (speed 2) (safety 0)))
|
|
804
|
813
|
(labels
|
| ... |
... |
@@ -982,13 +991,13 @@ |
|
982
|
991
|
(((complex double-double-float)
|
|
983
|
992
|
(foreach (complex rational) (complex single-float) (complex double-float)
|
|
984
|
993
|
(complex double-double-float)))
|
|
985
|
|
- (c::cdiv-double-double-float x
|
|
|
994
|
+ (cdiv-double-double-float x
|
|
986
|
995
|
(coerce y '(complex double-double-float))))
|
|
987
|
996
|
|
|
988
|
997
|
(((foreach integer ratio single-float double-float double-double-float
|
|
989
|
998
|
(complex rational) (complex single-float) (complex double-float))
|
|
990
|
999
|
(complex double-double-float))
|
|
991
|
|
- (c::cdiv-double-double-float (coerce x '(complex double-double-float))
|
|
|
1000
|
+ (cdiv-double-double-float (coerce x '(complex double-double-float))
|
|
992
|
1001
|
y))
|
|
993
|
1002
|
|
|
994
|
1003
|
(((foreach integer ratio single-float double-float double-double-float)
|