Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl

Commits:

3 changed files:

Changes:

  • src/code/numbers.lisp
    ... ... @@ -594,6 +594,7 @@
    594 594 
     	    (build-ratio (truncate x gcd) (truncate y gcd))))))
    595 595
    596 596
    597 
    +#+nil
    597 598 
     (defun two-arg-/ (x y)
    598 599 
       (number-dispatch ((x number) (y number))
    599 600 
         (float-contagion / x y (ratio integer))
    ... ... @@ -656,6 +657,441 @@
    656 657 
            (build-ratio (maybe-truncate nx gcd)
    657 658 
     		    (* (maybe-truncate y gcd) (denominator x)))))))
    658 659
    660 
    +;; An implementation of Baudin and Smith's robust complex division.
    661 
    +;; This is a pretty straightforward translation of the original in
    662 
    +;; https://arxiv.org/pdf/1210.4539.
    663 
    +#+nil
    664 
    +(eval-when (:compile-toplevel :load-toplevel :execute)
    665 
    +(macrolet
    666 
    +    ((frob (type max-float min-float eps-float)
    667 
    +       (let ((name (symbolicate "CDIV-" type)))
    668 
    +	 `(progn
    669 
    +	    (let ((ov ,max-float)
    670 
    +		  (un ,min-float)
    671 
    +		  (eps ,eps-float))
    672 
    +	      (defun ,name (x y)
    673 
    +		(declare (type (complex ,type) x y))
    674 
    +		(labels
    675 
    +		    ((internal-compreal (a b c d r tt)
    676 
    +		       (declare (double-float a b c d r tt))
    677 
    +		       ;; Iteration 1:  compare r against DBL_MIN instead of 0
    678 
    +		       (cond ((/= r 0)
    679 
    +			      (let ((br (* b r)))
    680 
    +				(if (/= br 0)
    681 
    +				    (* (+ a br) tt)
    682 
    +				    (+ (* a tt)
    683 
    +				       (* (* b tt)
    684 
    +					  r)))))
    685 
    +			     (t
    686 
    +			      (* (+ a (* d (/ b c)))
    687 
    +				 tt))))
    688 
    +
    689 
    +		     (robust-subinternal (a b c d)
    690 
    +		       (declare (double-float a b c d))
    691 
    +		       (let* ((r (/ d c))
    692 
    +			      (tt (/ (+ c (* d r)))))
    693 
    +			 (let ((e (internal-compreal a b c d r tt))
    694 
    +			       (f (internal-compreal b (- a) c d r tt)))
    695 
    +			   #+nil(format t "subint: e f = ~A ~A~%" e f)
    696 
    +			   (values e
    697 
    +				   f))))
    698 
    +
    699 
    +		     (robust-internal (x y)
    700 
    +		       (declare (type (complex ,type) x y))
    701 
    +		       #+nil(format t "x = ~A, y = ~A~%" x y)
    702 
    +		       (let ((a (realpart x))
    703 
    +			     (b (imagpart x))
    704 
    +			     (c (realpart y))
    705 
    +			     (d (imagpart y)))
    706 
    +			 (cond
    707 
    +			   ((<= (abs d) (abs c))
    708 
    +			    (multiple-value-bind (e f)
    709 
    +				(robust-subinternal a b c d)
    710 
    +			      #+nil(format t "1: e f = ~A ~A~%" e f)
    711 
    +			      (complex e f)))
    712 
    +			   (t
    713 
    +			    (multiple-value-bind (e f)
    714 
    +				(robust-subinternal b a d c)
    715 
    +			      #+nil(format t "2: e f = ~A ~A~%" e f)
    716 
    +			      (complex e (- f))))))))
    717 
    +		  (let* ((a (realpart x))
    718 
    +			 (b (imagpart x))
    719 
    +			 (c (realpart y))
    720 
    +			 (d (imagpart y))
    721 
    +			 (ab (max (abs a) (abs b)))
    722 
    +			 (cd (max (abs c) (abs d)))
    723 
    +			 (b 2d0)
    724 
    +			 (s 1d0)
    725 
    +			 (be (/ b (* eps eps))))
    726 
    +		    (when (>= ab (/ ov 2))
    727 
    +		      (setf x (/ x 2)
    728 
    +			    s (* s 2)))
    729 
    +		    (when (>= cd (/ ov 2))
    730 
    +		      (setf y (/ y 2)
    731 
    +			    s (/ s 2)))
    732 
    +		    (when (<= ab (* un (/ b eps)))
    733 
    +		      (setf x (* x be)
    734 
    +			    s (/ s be)))
    735 
    +		    (when (<= cd (* un (/ b eps)))
    736 
    +		      (setf y (* y be)
    737 
    +			    s (* s be)))
    738 
    +		    ;; Iteration 2
    739 
    +		    #+nil
    740 
    +		    (when (< cd eps)
    741 
    +		      (setf x (/ x eps))
    742 
    +		      (setf y (/ y eps)))
    743 
    +
    744 
    +		    ;; Iteration 4
    745 
    +		    #+nil
    746 
    +		    (when (> cd (/ ov 2))
    747 
    +		      (setf x (/ x 2)
    748 
    +			    y (/ y 2)))
    749 
    +
    750 
    +		    (* s
    751 
    +		       (robust-internal x y))))))))))
    752 
    +  ;; The eps value for double-float is determined by looking at the
    753 
    +  ;; value of %eps in Scilab.
    754 
    +  (frob double-float most-positive-double-float least-positive-normalized-double-float
    755 
    +	(scale-float 1d0 -52))
    756 
    +  ;; The eps value here is a guess.
    757 
    +  (frob single-float most-positive-single-float least-positive-normalized-single-float
    758 
    +	(scale-float 1f0 -22)))
    759 
    +)
    760 
    +
    761 
    +(let ((ov most-positive-double-float)
    762 
    +      (un least-positive-normalized-double-float)
    763 
    +      (eps (scale-float 1d0 -52)))
    764 
    +  (defun cdiv-double-float (x y)
    765 
    +    (declare (type (complex double-float) x y))
    766 
    +    (labels
    767 
    +	((internal-compreal (a b c d r tt)
    768 
    +	   (declare (double-float a b c d r tt))
    769 
    +	   ;; Iteration 1:  compare r against DBL_MIN instead of 0
    770 
    +	   (cond ((/= r 0)
    771 
    +		  (let ((br (* b r)))
    772 
    +		    (if (/= br 0)
    773 
    +			(* (+ a br) tt)
    774 
    +			(+ (* a tt)
    775 
    +			   (* (* b tt)
    776 
    +			      r)))))
    777 
    +		 (t
    778 
    +		  (* (+ a (* d (/ b c)))
    779 
    +		     tt))))
    780 
    +
    781 
    +	 (robust-subinternal (a b c d)
    782 
    +	   (declare (double-float a b c d))
    783 
    +	   (let* ((r (/ d c))
    784 
    +		  (tt (/ (+ c (* d r)))))
    785 
    +	     (let ((e (internal-compreal a b c d r tt))
    786 
    +		   (f (internal-compreal b (- a) c d r tt)))
    787 
    +	       #+nil(format t "subint: e f = ~A ~A~%" e f)
    788 
    +	       (values e
    789 
    +		       f))))
    790 
    +
    791 
    +	 (robust-internal (x y)
    792 
    +	   (declare (type (complex double-float) x y))
    793 
    +	   #+nil(format t "x = ~A, y = ~A~%" x y)
    794 
    +	   (let ((a (realpart x))
    795 
    +		 (b (imagpart x))
    796 
    +		 (c (realpart y))
    797 
    +		 (d (imagpart y)))
    798 
    +	     (cond
    799 
    +	       ((<= (abs d) (abs c))
    800 
    +		(multiple-value-bind (e f)
    801 
    +		    (robust-subinternal a b c d)
    802 
    +		  #+nil(format t "1: e f = ~A ~A~%" e f)
    803 
    +		  (complex e f)))
    804 
    +	       (t
    805 
    +		(multiple-value-bind (e f)
    806 
    +		    (robust-subinternal b a d c)
    807 
    +		  #+nil(format t "2: e f = ~A ~A~%" e f)
    808 
    +		  (complex e (- f))))))))
    809 
    +      (let* ((a (realpart x))
    810 
    +	     (b (imagpart x))
    811 
    +	     (c (realpart y))
    812 
    +	     (d (imagpart y))
    813 
    +	     (ab (max (abs a) (abs b)))
    814 
    +	     (cd (max (abs c) (abs d)))
    815 
    +	     (b 2d0)
    816 
    +	     (s 1d0)
    817 
    +	     (be (/ b (* eps eps))))
    818 
    +	(when (>= ab (/ ov 2))
    819 
    +	  (setf x (/ x 2)
    820 
    +		s (* s 2)))
    821 
    +	(when (>= cd (/ ov 2))
    822 
    +	  (setf y (/ y 2)
    823 
    +		s (/ s 2)))
    824 
    +	(when (<= ab (* un (/ b eps)))
    825 
    +	  (setf x (* x be)
    826 
    +		s (/ s be)))
    827 
    +	(when (<= cd (* un (/ b eps)))
    828 
    +	  (setf y (* y be)
    829 
    +		s (* s be)))
    830 
    +	;; Iteration 2
    831 
    +	#+nil
    832 
    +	(when (< cd eps)
    833 
    +	  (setf x (/ x eps))
    834 
    +	  (setf y (/ y eps)))
    835 
    +
    836 
    +	;; Iteration 4
    837 
    +	#+nil
    838 
    +	(when (> cd (/ ov 2))
    839 
    +	  (setf x (/ x 2)
    840 
    +		y (/ y 2)))
    841 
    +
    842 
    +	(* s
    843 
    +	   (robust-internal x y))))))
    844 
    +
    845 
    +(defun cdiv-single-float (x y)
    846 
    +  (declare (type (complex single-float) x y))
    847 
    +  (let ((a (float (realpart x) 1d0))
    848 
    +	(b (float (imagpart x) 1d0))
    849 
    +	(c (float (realpart y) 1d0))
    850 
    +	(d (float (imagpart y) 1d0)))
    851 
    +    (cond ((< (abs c) (abs d))
    852 
    +	   (let* ((r (/ c d))
    853 
    +		  (denom (+ (* c r) d))
    854 
    +		  (e (float (/ (+ (* a r) b) denom) 1f0))
    855 
    +		  (f (float (/ (- (* b r) a) denom) 1f0)))
    856 
    +	     (complex e f)))
    857 
    +	  (t
    858 
    +	   (let* ((r (/ d c))
    859 
    +		  (denom (+ c (* d r)))
    860 
    +		  (e (float (/ (+ a (* b r)) denom) 1f0))
    861 
    +		  (f (float (/ (- b (* a r)) denom) 1f0)))
    862 
    +	     (complex e f))))))
    863 
    +
    864 
    +(defun two-arg-/ (x y)
    865 
    +  (number-dispatch ((x number) (y number))
    866 
    +    (float-contagion / x y (ratio integer))
    867 
    +
    868 
    +    (((complex single-float)
    869 
    +      (foreach (complex rational) (complex single-float)))
    870 
    +     (cdiv-single-float x (coerce y '(complex single-float))))
    871 
    +    (((complex double-float)
    872 
    +      (foreach (complex rational) (complex single-float) (complex double-float)))
    873 
    +     (cdiv-double-float x (coerce y '(complex double-float))))
    874 
    +
    875 
    +    (((foreach integer ratio single-float (complex rational))
    876 
    +      (complex single-float))
    877 
    +     (cdiv-single-float (coerce x '(complex single-float))
    878 
    +			y))
    879 
    +
    880 
    +    (((foreach integer ratio single-float double-float (complex rational) (complex single-float))
    881 
    +      (complex double-float))
    882 
    +     (cdiv-double-float (coerce x '(complex double-float))
    883 
    +			y))
    884 
    +    (((complex double-double-float)
    885 
    +      (foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float)))
    886 
    +     (let ((a (realpart x))
    887 
    +	   (b (imagpart x))
    888 
    +	   (c (realpart y))
    889 
    +	   (d (imagpart y)))
    890 
    +       (cond ((< (abs c) (abs d))
    891 
    +	      (let* ((r (/ c d))
    892 
    +		     (denom (+ (* c r) d))
    893 
    +		     (e (/ (+ (* a r) b) denom))
    894 
    +		     (f (/ (- (* b r) a) denom)))
    895 
    +		(canonical-complex e f)))
    896 
    +	     (t
    897 
    +	      (let* ((r (/ d c))
    898 
    +		     (denom (+ c (* d r)))
    899 
    +		     (e (/ (+ a (* b r)) denom))
    900 
    +		     (f (/ (- b (* a r)) denom)))
    901 
    +		(canonical-complex e f))))))
    902 
    +
    903 
    +    (((foreach integer ratio single-float double-float double-double-float
    904 
    +	       (complex rational) (complex single-float) (complex double-float))
    905 
    +      (complex double-double-float))
    906 
    +     (let ((a (realpart x))
    907 
    +	   (b (imagpart x))
    908 
    +	   (c (realpart y))
    909 
    +	   (d (imagpart y)))
    910 
    +       (cond ((< (abs c) (abs d))
    911 
    +	      (let* ((r (/ c d))
    912 
    +		     (denom (+ (* c r) d))
    913 
    +		     (e (/ (+ (* a r) b) denom))
    914 
    +		     (f (/ (- (* b r) a) denom)))
    915 
    +		(canonical-complex e f)))
    916 
    +	     (t
    917 
    +	      (let* ((r (/ d c))
    918 
    +		     (denom (+ c (* d r)))
    919 
    +		     (e (/ (+ a (* b r)) denom))
    920 
    +		     (f (/ (- b (* a r)) denom)))
    921 
    +		(canonical-complex e f))))))
    922 
    +
    923 
    +    (((foreach integer ratio single-float double-float double-double-float)
    924 
    +      (complex rational))
    925 
    +     (let* ((ry (realpart y))
    926 
    +	    (iy (imagpart y)))
    927 
    +       (if (> (abs ry) (abs iy))
    928 
    +	   (let* ((r (/ iy ry))
    929 
    +		  (dn (* ry (+ 1 (* r r)))))
    930 
    +	     (canonical-complex (/ x dn)
    931 
    +				(/ (- (* x r)) dn)))
    932 
    +	   (let* ((r (/ ry iy))
    933 
    +		  (dn (* iy (+ 1 (* r r)))))
    934 
    +	     (canonical-complex (/ (* x r) dn)
    935 
    +				(/ (- x) dn))))))
    936 
    +    (((complex rational)
    937 
    +      (complex rational))
    938 
    +     (let ((a (realpart x))
    939 
    +	   (b (imagpart x))
    940 
    +	   (c (realpart y))
    941 
    +	   (d (imagpart y)))
    942 
    +       (cond ((< (abs c) (abs d))
    943 
    +	      (let* ((r (/ c d))
    944 
    +		     (denom (+ (* c r) d))
    945 
    +		     (e (/ (+ (* a r) b) denom))
    946 
    +		     (f (/ (- (* b r) a) denom)))
    947 
    +		(canonical-complex e f)))
    948 
    +	     (t
    949 
    +	      (let* ((r (/ d c))
    950 
    +		     (denom (+ c (* d r)))
    951 
    +		     (e (/ (+ a (* b r)) denom))
    952 
    +		     (f (/ (- b (* a r)) denom)))
    953 
    +		(canonical-complex e f))))))
    954 
    +
    955 
    +    (((foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float))
    956 
    +      (or rational float))
    957 
    +     (canonical-complex (/ (realpart x) y)
    958 
    +			(/ (imagpart x) y)))
    959 
    +
    960 
    +    ((double-float
    961 
    +      (complex single-float))
    962 
    +     (cdiv-double-float (coerce x '(complex double-float))
    963 
    +			(coerce y '(complex double-float))))
    964 
    +
    965 
    +    #+nil
    966 
    +    (((foreach integer ratio single-float double-float) complex)
    967 
    +     (let* ((ry (realpart y))
    968 
    +	    (iy (imagpart y)))
    969 
    +       (if (> (abs ry) (abs iy))
    970 
    +	   (let* ((r (/ iy ry))
    971 
    +		  (dn (* ry (+ 1 (* r r)))))
    972 
    +	     (canonical-complex (/ x dn)
    973 
    +				(/ (- (* x r)) dn)))
    974 
    +	   (let* ((r (/ ry iy))
    975 
    +		  (dn (* iy (+ 1 (* r r)))))
    976 
    +	     (canonical-complex (/ (* x r) dn)
    977 
    +				(/ (- x) dn))))))
    978 
    +    #+nil
    979 
    +    ((complex (or rational float))
    980 
    +     (canonical-complex (/ (realpart x) y)
    981 
    +			(/ (imagpart x) y)))
    982 
    +
    983 
    +    ((ratio ratio)
    984 
    +     (let* ((nx (numerator x))
    985 
    +	    (dx (denominator x))
    986 
    +	    (ny (numerator y))
    987 
    +	    (dy (denominator y))
    988 
    +	    (g1 (gcd nx ny))
    989 
    +	    (g2 (gcd dx dy)))
    990 
    +       (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
    991 
    +		    (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))
    992 
    +
    993 
    +    ((integer integer)
    994 
    +     (integer-/-integer x y))
    995 
    +
    996 
    +    ((integer ratio)
    997 
    +     (if (zerop x)
    998 
    +	 0
    999 
    +	 (let* ((ny (numerator y))
    1000 
    +		(dy (denominator y))
    1001 
    +		(gcd (gcd x ny)))
    1002 
    +	   (build-ratio (* (maybe-truncate x gcd) dy)
    1003 
    +			(maybe-truncate ny gcd)))))
    1004 
    +
    1005 
    +    ((ratio integer)
    1006 
    +     (let* ((nx (numerator x))
    1007 
    +	    (gcd (gcd nx y)))
    1008 
    +       (build-ratio (maybe-truncate nx gcd)
    1009 
    +		    (* (maybe-truncate y gcd) (denominator x)))))))
    1010 
    +
    1011 
    +
    1012 
    +#+nil
    1013 
    +(defun two-arg-/ (x y)
    1014 
    +  (number-dispatch ((x number) (y number))
    1015 
    +    (float-contagion / x y (ratio integer))
    1016 
    +
    1017 
    +    (((complex double-float) (complex double-float))
    1018 
    +     (cdiv-double-float x y))
    1019 
    +    (((complex double-float) (complex single-float))
    1020 
    +     (cdiv-double-float x (coerce y '(complex double-float))))
    1021 
    +
    1022 
    +    (((complex single-float) (complex double-float))
    1023 
    +     (cdiv-double-float (coerce y '(complex double-float))
    1024 
    +			y))
    1025 
    +    (((complex single-float) (complex single-float))
    1026 
    +     (cdiv-single-float x y))
    1027 
    +
    1028 
    +    #+nil
    1029 
    +    (((foreach (complex single-float) (complex double-float))
    1030 
    +      (complex double-float))
    1031 
    +     (cdiv-double-float (complex (float (realpart x) 1d0)
    1032 
    +				 (float (imagpart x) 1d0))
    1033 
    +			y))
    1034 
    +    (((foreach (complex rational) (complex double-double-float))
    1035 
    +      (foreach (complex rational) (complex double-double-float)))
    1036 
    +     (let* ((rx (realpart x))
    1037 
    +	    (ix (imagpart x))
    1038 
    +	    (ry (realpart y))
    1039 
    +	    (iy (imagpart y)))
    1040 
    +       (if (> (abs ry) (abs iy))
    1041 
    +	   (let* ((r (/ iy ry))
    1042 
    +		  (dn (+ ry (* r iy))))
    1043 
    +	     (canonical-complex (/ (+ rx (* ix r)) dn)
    1044 
    +				(/ (- ix (* rx r)) dn)))
    1045 
    +	   (let* ((r (/ ry iy))
    1046 
    +		  (dn (+ iy (* r ry))))
    1047 
    +	     (canonical-complex (/ (+ (* rx r) ix) dn)
    1048 
    +				(/ (- (* ix r) rx) dn))))))
    1049 
    +    (((foreach integer ratio single-float double-float)
    1050 
    +      (foreach (complex rational) (complex double-double-float)))
    1051 
    +     (let* ((ry (realpart y))
    1052 
    +	    (iy (imagpart y)))
    1053 
    +       (if (> (abs ry) (abs iy))
    1054 
    +	   (let* ((r (/ iy ry))
    1055 
    +		  (dn (* ry (+ 1 (* r r)))))
    1056 
    +	     (canonical-complex (/ x dn)
    1057 
    +				(/ (- (* x r)) dn)))
    1058 
    +	   (let* ((r (/ ry iy))
    1059 
    +		  (dn (* iy (+ 1 (* r r)))))
    1060 
    +	     (canonical-complex (/ (* x r) dn)
    1061 
    +				(/ (- x) dn))))))
    1062 
    +    (((foreach (complex rational) (complex double-double-float))
    1063 
    +      (or rational float))
    1064 
    +     (canonical-complex (/ (realpart x) y)
    1065 
    +			(/ (imagpart x) y)))
    1066 
    +
    1067 
    +    ((ratio ratio)
    1068 
    +     (let* ((nx (numerator x))
    1069 
    +	    (dx (denominator x))
    1070 
    +	    (ny (numerator y))
    1071 
    +	    (dy (denominator y))
    1072 
    +	    (g1 (gcd nx ny))
    1073 
    +	    (g2 (gcd dx dy)))
    1074 
    +       (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
    1075 
    +		    (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))
    1076 
    +
    1077 
    +    ((integer integer)
    1078 
    +     (integer-/-integer x y))
    1079 
    +
    1080 
    +    ((integer ratio)
    1081 
    +     (if (zerop x)
    1082 
    +	 0
    1083 
    +	 (let* ((ny (numerator y))
    1084 
    +		(dy (denominator y))
    1085 
    +		(gcd (gcd x ny)))
    1086 
    +	   (build-ratio (* (maybe-truncate x gcd) dy)
    1087 
    +			(maybe-truncate ny gcd)))))
    1088 
    +
    1089 
    +    ((ratio integer)
    1090 
    +     (let* ((nx (numerator x))
    1091 
    +	    (gcd (gcd nx y)))
    1092 
    +       (build-ratio (maybe-truncate nx gcd)
    1093 
    +		    (* (maybe-truncate y gcd) (denominator x)))))))
    1094 
    +
    659 1095
    660 1096 
     (defun %negate (n)
    661 1097 
       (number-dispatch ((n number))

  • src/compiler/float-tran.lisp
    ... ... @@ -1917,104 +1917,13 @@
    1917 1917 
       (frob (complex single-float) 1f0)
    1918 1918 
       (frob (complex double-float) 1d0))
    1919 1919
    1920 
    -;; An implementation of Baudin and Smith's robust complex division.
    1921 
    -;; This is a pretty straightforward translation of the original in
    1922 
    -;; https://arxiv.org/pdf/1210.4539.
    1923 1920 
     (macrolet
    1924 
    -    ((frob (type max-float min-float eps-float)
    1925 
    -       (let ((name (symbolicate "CDIV-" type)))
    1926 
    -	 `(progn
    1927 
    -	    (let ((ov ,max-float)
    1928 
    -		  (un ,min-float)
    1929 
    -		  (eps ,eps-float))
    1930 
    -	      (defun ,name (x y)
    1931 
    -		(declare (type (complex ,type) x y))
    1932 
    -		(labels
    1933 
    -		    ((internal-compreal (a b c d r tt)
    1934 
    -		       (declare (,type a b c d r tt))
    1935 
    -		       (cond ((zerop r)
    1936 
    -			      (let ((br (* b r)))
    1937 
    -				(if (/= br 0)
    1938 
    -				    (* (+ a br) tt)
    1939 
    -				    (+ (* a tt)
    1940 
    -				       (* (* b tt)
    1941 
    -					  r)))))
    1942 
    -			     (t
    1943 
    -			      (* (+ a (* d (/ b c)))
    1944 
    -				 tt))))
    1945 
    -
    1946 
    -		     (robust-subinternal (a b c d)
    1947 
    -		       (declare (,type a b c d))
    1948 
    -		       (let* ((r (/ d c))
    1949 
    -			      (tt (/ (+ c (* d r)))))
    1950 
    -			 (values (internal-compreal a b c d r tt)
    1951 
    -				 (internal-compreal b (- a) c d r tt))))
    1952 
    -
    1953 
    -		     (robust-internal (x y)
    1954 
    -		       (declare (type (complex ,type) x y))
    1955 
    -		       (let ((a (realpart x))
    1956 
    -			     (b (imagpart x))
    1957 
    -			     (c (realpart y))
    1958 
    -			     (d (imagpart y)))
    1959 
    -			 (cond
    1960 
    -			   ((<= (abs d) (abs c))
    1961 
    -			    (multiple-value-bind (e f)
    1962 
    -				(robust-subinternal a b c d)
    1963 
    -			      (complex e f)))
    1964 
    -			   (t
    1965 
    -			    (multiple-value-bind (e f)
    1966 
    -				(robust-subinternal b a d c)
    1967 
    -			      (complex e (- f))))))))
    1968 
    -		  (let* ((a (realpart x))
    1969 
    -			 (b (imagpart x))
    1970 
    -			 (c (realpart y))
    1971 
    -			 (d (imagpart y))
    1972 
    -			 (ab (max (abs a) (abs b)))
    1973 
    -			 (cd (max (abs c) (abs d)))
    1974 
    -			 (b 2d0)
    1975 
    -			 (s 1d0)
    1976 
    -			 (be (/ b (* eps eps))))
    1977 
    -		    (when (>= ab (/ ov 2))
    1978 
    -		      (setf x (/ x 2)
    1979 
    -			    s (* s 2)))
    1980 
    -		    (when (>= cd (/ ov 2))
    1981 
    -		      (setf y (/ y 2)
    1982 
    -			    s (/ s 2)))
    1983 
    -		    (when (<= ab (* un (/ b eps)))
    1984 
    -		      (setf x (* x be)
    1985 
    -			    s (/ s be)))
    1986 
    -		    (when (<= cd (* un (/ b eps)))
    1987 
    -		      (setf y (* y be)
    1988 
    -			    s (* s be)))
    1989 
    -		    ;; Iteration 2
    1990 
    -		    #+nil
    1991 
    -		    (when (< cd eps)
    1992 
    -		      (setf x (/ x eps))
    1993 
    -		      (setf y (/ y eps)))
    1994 
    -
    1995 
    -		    ;; Iteration 4
    1996 
    -		    #+nil
    1997 
    -		    (when (> cd (/ ov 2))
    1998 
    -		      (setf x (/ x 2)
    1999 
    -			    y (/ y 2)))
    2000 
    -
    2001 
    -		    (* s
    2002 
    -		       (robust-internal x y))))))))))
    2003 
    -  ;; The eps value for double-float is determined by looking at the
    2004 
    -  ;; value of %eps in Scilab.
    2005 
    -  (frob double-float most-positive-double-float least-positive-normalized-double-float
    2006 
    -	(scale-float 1d0 -52))
    2007 
    -  ;; The eps value here is a guess.
    2008 
    -  (frob single-float most-positive-single-float least-positive-normalized-single-float
    2009 
    -	(scale-float 1f0 -22)))
    2010 
    -
    2011 
    -(macrolet
    2012 
    -    ((frob (type)
    2013 
    -       (let ((name (symbolicate "CDIV-" type)))
    1921 
    +    ((frob (type name)
    2014 1922 
            `(deftransform / ((x y) ((complex ,type) (complex ,type)) *)
    2015 
    -		      (,name x y)))))
    2016 
    -  (frob double-float)
    2017 
    -  (frob single-float))
    1923 
    +		      (,name x y))))
    1924 
    +  (frob double-float kernel::cdiv-double-float)
    1925 
    +  (frob single-float kernel::cdiv-single-float))
    1926 
    +
    2018 1927
    2019 1928 
     ;;;; Complex contagion:
    2020 1929

  • tests/float.lisp
    ... ... @@ -388,7 +388,7 @@
    388 388 
        (list (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074))
    389 389 
     	 (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074))
    390 390 
     	 (complex 0.6d0 0.2d0)
    391 
    -	 53 least-positive-double-float)
    391 
    +	 52 least-positive-double-float)
    392 392 
        ;; 9
    393 393 
        (list (complex (scale-float 1d0 1015) (scale-float 1d0 -989))
    394 394 
     	 (complex (scale-float 1d0 1023) (scale-float 1d0 1023))
    ... ... @@ -456,7 +456,7 @@
    456 456 
           (min (real-ulp (realpart diff))
    457 457 
     	   (real-ulp (imagpart diff))))))
    458 458
    459 
    -(define-test complex-division
    459 
    +(define-test complex-division.double
    460 460 
       (:tag :issues)
    461 461 
       (loop for k from 1
    462 462 
     	for test in *test-cases*
    ... ... @@ -473,3 +473,35 @@
    473 473 
     			     k x y z z-true diff rel)
    474 474 
     	       (assert-true (<= ulp max-ulp)
    475 475 
     			    k x y z z-true diff ulp max-ulp)))))
    476 
    +
    477 
    +(define-test complex-division.misc
    478 
    +    (:tag :issue)
    479 
    +  (let ((num '(1
    480 
    +	       1/2
    481 
    +	       1.0
    482 
    +	       1d0
    483 
    +	       #c(1 2)
    484 
    +	       #c(1.0 2.0)
    485 
    +	       #c(1d0 2d0)
    486 
    +	       #c(1w0 2w0))))
    487 
    +    ;; Try all combinations of divisions of different types.  This is
    488 
    +    ;; primarily to test that we got all the numeric contagion cases
    489 
    +    ;; for division in CL:/.
    490 
    +    (dolist (x num)
    491 
    +      (dolist (y num)
    492 
    +	(assert-true (/ x y)
    493 
    +		     x y)))))
    494 
    +
    495 
    +(define-test complex-division.single
    496 
    +    (:tag :issues)
    497 
    +  (let ((x #c(1 2))
    498 
    +	(y (complex (expt 2 127) (expt 2 127)))
    499 
    +	(expected (coerce (/ x y)
    500 
    +			  '(complex single-float))))
    501 
    +    ;; A naive implementation of complex division would cause an
    502 
    +    ;; overflow in computing the denominator.
    503 
    +    (assert-equal expected
    504 
    +		  (/ (coerce x '(complex single-float))
    505 
    +		     (coerce y '(complex single-float)))
    506 
    +		  x
    507 
    +		  y)))

    • —
    View it on GitLab.
    You're receiving this email because of your account on gitlab.common-lisp.net. Manage all notifications · Help Notification message regarding https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/e11ed3eb2b57116cf92220cd470829eb2e98143d...825f7d96ec6bc182a2e6a557702818220211e7f6 at 1765408966