Changeset 15627
 Timestamp:
 Jan 31, 2013, 5:03:29 PM (6 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/source/level0/l0float.lisp
r15553 r15627 23 23 (require "NUMBERMACROS") 24 24 (require :numbercasemacro) 25 (defconstant singlefloatpi (coerce pi 'singlefloat)) 26 (defconstant doublefloathalfpi (/ pi 2)) 27 (defconstant singlefloathalfpi (coerce (/ pi 2) 'singlefloat)) 28 (defconstant singlefloatlog2 0.6931472) ; (log 2) 29 (defconstant doublefloatlog2 0.6931471805599453d0) ; (log 2.0d0) 30 (defconstant doublefloatlog2^23 15.942385152878742d0) ; (log (expt 2 23)) 25 31 ) 26 32 … … 729 735 730 736 731 737 ;;; Multiply by i (with additional optional scale factor) 738 ;;; Does the "right thing" with minus zeroes (see CLTL2) 739 (defun i* (number &optional (scale 1)) 740 (complex (* ( scale) (imagpart number)) 741 (* scale (realpart number)))) 742 743 ;;; complex atanh 744 (defun %complexatanh (z) 745 (let* ((rx (realpart z)) 746 (ix (imagpart z)) 747 (sign (typecase rx 748 (doublefloat (%doublefloatsign rx)) 749 (shortfloat (%shortfloatsign rx)) 750 (t (minusp rx)))) 751 (x rx) 752 (y ix) 753 (y1 (abs y)) 754 ra ia) 755 ;; following code requires nonnegative x 756 (when sign 757 (setf x ( x)) 758 (setf y ( y))) 759 (cond ((> (max x y1) 1.8014399e+16) 760 ;; large value escape 761 (setq ra (if (> x y1) 762 (let ((r (/ y x))) 763 (/ (/ x) (1+ (* r r)))) 764 (let ((r (/ x y))) 765 (/ (/ r y) (1+ (* r r)))))) 766 (setq ia (typecase y 767 (doublefloat (floatsign y doublefloathalfpi)) 768 (singlefloat (floatsign y singlefloathalfpi)) 769 (t (if (minusp y) #.( singlefloathalfpi) singlefloathalfpi))))) 770 ((= x 1) 771 (setq ra (if (< y1 1e9) 772 (/ (loge (/ 2 y1)) 2) 773 (/ (log1+ (/ 4 (* y y))) 4))) 774 (setq ia (/ (atan (/ 2 y) 1) 2))) 775 (t 776 (let ((r2 (+ (* x x) (* y y)))) 777 (setq ra (/ (log1+ (/ (* 4 x) (1+ (+ (* 2 x) r2)))) 4)) 778 (setq ia (/ (atan (* 2 y) ( 1 r2)) 2))))) 779 ;; fixup signs, with special case for real arguments 780 (cond (sign 781 (setq ra ( ra)) 782 (when (typep z 'complex) 783 (setq ia ( ia)))) 784 (t 785 (unless (typep z 'complex) 786 (setq ia ( ia))))) 787 (complex ra ia))) 732 788 733 789 (defun atan (y &optional (x nil xp)) 734 790 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied." 735 (if xp 736 (if (or (typep x 'doublefloat) 737 (typep y 'doublefloat)) 738 (withstackdoublefloats ((dy y) 739 (dx x)) 740 (%dfatan2 dy dx)) 741 #+32bittarget 742 (target::withstackshortfloats ((sy y) 743 (sx x)) 744 (%sfatan2! sy sx)) 745 #+64bittarget 746 (%sfatan2 (%shortfloat y) (%shortfloat x))) 747 (if (typep y 'complex) 748 (let ((iy (complex ( (imagpart y)) (realpart y)))) 749 (/ ( (log (+ 1 iy)) (log ( 1 iy))) 750 #c(0 2))) 751 (if (typep y 'doublefloat) 752 (%doublefloatatan! y (%makedfloat)) 753 #+32bittarget 754 (target::withstackshortfloats ((sy y)) 755 (%singlefloatatan! sy (%makesfloat))) 756 #+64bittarget 757 (%singlefloatatan (%shortfloat y)) 758 )))) 791 (cond (xp 792 (cond ((or (typep x 'doublefloat) 793 (typep y 'doublefloat)) 794 (withstackdoublefloats ((dy y) 795 (dx x)) 796 (%dfatan2 dy dx))) 797 (t 798 (when (and (rationalp x) (rationalp y)) 799 ;; rescale arguments so that the maximum absolute value is 1 800 (let ((x1 (abs x)) (y1 (abs y))) 801 (cond ((> y1 x1) 802 (setf x (/ x y1)) 803 (setf y (signum y))) 804 ((not (zerop x)) 805 (setf y (/ y x1)) 806 (setf x (signum x)))))) 807 #+32bittarget 808 (target::withstackshortfloats ((sy y) 809 (sx x)) 810 (%sfatan2! sy sx)) 811 #+64bittarget 812 (%sfatan2 (%shortfloat y) (%shortfloat x))))) 813 ((typep y 'doublefloat) 814 (%doublefloatatan! y (%makedfloat))) 815 ((typep y 'singlefloat) 816 #+32bittarget 817 (%singlefloatatan! y (%makesfloat)) 818 #+64bittarget 819 (%singlefloatatan y)) 820 ((typep y 'rational) 821 (cond ((<= (abs y) mostpositiveshortfloat) 822 #+32bittarget 823 (target::withstackshortfloats ((sy y)) 824 (%singlefloatatan! sy (%makesfloat))) 825 #+64bittarget 826 (%singlefloatatan (%shortfloat y))) 827 ((minusp y) 828 #.( singlefloathalfpi)) 829 (t 830 singlefloathalfpi))) 831 (t 832 (let ((r (realpart y)) 833 (i (imagpart y))) 834 (if (zerop i) 835 (complex (atan r) i) 836 (i* (%complexatanh (complex ( i) r)) 1)))))) 759 837 760 838 … … 766 844 (if (zerop x) 767 845 (reportbadarg x '(not (satisfies zerop) )) 768 (if (floatp x) (float 0.0d0 x) 0)) 769 (/ (loge x) (loge b))) 846 ;; ** CORRECT THE CONTAGION for complex args ** 847 (+ 0 (* x b))) 848 ;; do the float/rational contagion before the division 849 ;; but take care with negative zeroes 850 (let ((x1 (realpart x)) 851 (b1 (realpart b))) 852 (if (and (typep x1 'float) 853 (typep b1 'float)) 854 (/ (loge (* x (float 1.0 b1))) 855 (loge (* b (float 1.0 x1)))) 856 (let ((r (/ (cond ((typep x 'rational) 857 (%rationallog x 1.0d0)) 858 ((typep x1 'rational) 859 (%rationalcomplexlog x 1.0d0)) 860 (t 861 (loge (* x 1.0d0)))) 862 (cond ((typep b 'rational) 863 (%rationallog b 1.0d0)) 864 ((typep b1 'rational) 865 (%rationalcomplexlog b 1.0d0)) 866 (t 867 (loge (* b 1.0d0))))))) 868 (cond ((or (typep x1 'doublefloat) 869 (typep b1 'doublefloat)) 870 r) 871 ((complexp r) 872 (complex (%shortfloat (realpart r)) 873 (%shortfloat (imagpart r)))) 874 (t 875 (%shortfloat r))))))) 770 876 (loge x))) 771 877 878 879 772 880 (defun loge (x) 773 (cond 774 ((bignump x) 775 (if (minusp x) 776 (complex (loge ( x)) pi) 777 (let* ((base1 3) 778 (guess (floor (1 (integerlength x)) 779 (log base1 2))) 780 (guess1 (* guess (loge base1)))) 781 (+ guess1 (loge (/ x (expt base1 guess))))))) 782 ((and (ratiop x) 783 (or (> x mostpositiveshortfloat) 784 (< x mostnegativeshortfloat))) 785 ( (loge (%numerator x)) (loge (%denominator x)))) 881 (cond 882 ((typep x 'doublefloat) 883 (if (%doublefloatsign x) 884 (withstackdoublefloats ((dx x)) 885 (complex (%doublefloatlog! (%%doublefloatabs! dx dx) (%makedfloat)) pi)) 886 (%doublefloatlog! x (%makedfloat)))) 887 ((typep x 'shortfloat) 888 #+32bittarget 889 (if (%shortfloatsign x) 890 (target::withstackshortfloats ((sx x)) 891 (complex (%singlefloatlog! (%%shortfloatabs! sx sx) (%makesfloat)) 892 singlefloatpi)) 893 (%singlefloatlog! x (%makesfloat))) 894 #+64bittarget 895 (if (%shortfloatsign x) 896 (complex (%singlefloatlog (%shortfloatabs (%shortfloat x))) 897 singlefloatpi) 898 (%singlefloatlog (%shortfloat x)))) 786 899 ((typep x 'complex) 787 (complex (loge (abs x)) (phase x))) 788 ((typep x 'doublefloat) 789 (withstackdoublefloats ((dx x)) 790 (if (minusp x) 791 (complex (%doublefloatlog! (%%doublefloatabs! dx dx) (%makedfloat)) pi) 792 (%doublefloatlog! dx (%makedfloat))))) 900 (if (typep (realpart x) 'rational) 901 (%rationalcomplexlog x 1.0s0) 902 ;; take care that intermediate results do not overflow/underflow: 903 ;; prescale argument by an appropriate power of two; 904 ;; we only need to scale for very large and very small values  905 ;; hence the various 'magic' numbers (values may not be too 906 ;; critical but do depend on the sqrt update to fix abs's operation) 907 (let ((m (max (abs (realpart x)) 908 (abs (imagpart x)))) 909 (logs 0) 910 (s 1)) 911 (if (typep m 'shortfloat) 912 (let ((expon ( (%shortfloatexp m) IEEEsinglefloatbias))) 913 (cond ((> expon 126) 914 (setq logs doublefloatlog2^23) 915 (setq s #.(ash 1 23))) 916 ((< expon 124) 917 (setq logs #.( doublefloatlog2^23)) 918 (setq s #.(/ 1.0s0 (ash 1 23)))))) 919 (let ((expon ( (%doublefloatexp m) IEEEdoublefloatbias))) 920 (cond ((> expon 1022) 921 (setq logs doublefloatlog2^23) 922 (setq s #.(ash 1 23))) 923 ((< expon 1020) 924 (setq logs #.( doublefloatlog2^23)) 925 (setq s #.(/ 1.0d0 (ash 1 23))))))) 926 (if (eql s 1) 927 (complex (logabs x) (phase x)) 928 (let ((temp (logabs (/ x s)))) 929 (complex (float (+ logs temp) temp) (phase x))))))) 793 930 (t 794 #+32bittarget 795 (target::withstackshortfloats ((sx x)) 796 (if (minusp x) 797 (complex (%singlefloatlog! (%%shortfloatabs! sx sx) (%makesfloat)) 798 #.(coerce pi 'shortfloat)) 799 (%singlefloatlog! sx (%makesfloat)))) 800 #+64bittarget 801 (if (minusp x) 802 (complex (%singlefloatlog (%shortfloatabs (%shortfloat x))) #.(coerce pi 'singlefloat)) 803 (%singlefloatlog (%shortfloat x))) 804 ))) 805 806 931 (%rationallog x 1.0s0)))) 932 933 ;;; helpers for rational log 934 (defun %rationallog (x prototype) 935 (cond ((minusp x) 936 (complex (%rationallog ( x) prototype) 937 (if (typep prototype 'shortfloat) 938 singlefloatpi 939 pi))) 940 ((bignump x) 941 ;(let* ((base1 3) 942 ; (guess (floor (1 (integerlength x)) 943 ; (log base1 2))) 944 ; (guess1 (* guess (loge base1)))) 945 ; (+ guess1 (loge (/ x (expt base1 guess))))) 946 ; Using base1 = 2 is *much* faster. Was there a reason for 3? 947 (let* ((guess (1 (integerlength x))) 948 (guess1 (* guess doublefloatlog2))) 949 (float (+ guess1 (loge (float (/ x (ash 1 guess)) 1.0d0))) prototype))) 950 ((and (ratiop x) 951 ;; Rational arguments near +1 can be specified with great precision: don't lose it 952 (cond ((< 0.5 x 3) 953 (log1+ (float ( x 1) prototype))) 954 ( 955 ;; Escape out small values as well as large 956 (or (> x mostpositiveshortfloat) 957 (< x leastpositivenormalizedshortfloat)) 958 ;; avoid loss of precision due to subtracting logs of numerator and denominator 959 (let* ((n (%numerator x)) 960 (d (%denominator x)) 961 (sn (1 (integerlength n))) 962 (sd (1 (integerlength d)))) 963 (float (+ (* ( sn sd) doublefloatlog2) 964 ( (log1+ (float (1 (/ n (ash 1 sn))) 1.0d0)) 965 (log1+ (float (1 (/ d (ash 1 sd))) 1.0d0)))) 966 prototype)))))) 967 ((typep prototype 'shortfloat) 968 #+32bittarget 969 (target::withstackshortfloats ((sx x)) 970 (%singlefloatlog! sx (%makesfloat))) 971 #+64bittarget 972 (%singlefloatlog (%shortfloat x))) 973 (t 974 (withstackdoublefloats ((dx x)) 975 (%doublefloatlog! dx (%makedfloat)))))) 976 977 ;;; (log1+ x) = (log (1+ x)) 978 ;;; but is much more precise near x = 0 979 (defun log1+ (x) 980 ;;(cond ((typep x 'complex) 981 ;; (let ((r (realpart x)) 982 ;; (i (imagpart x))) 983 ;; (if (and (< (abs r) 0.5) 984 ;; (< (abs i) 3)) 985 ;; (let* ((n (+ (* r (+ 2 r)) (* i i))) 986 ;; (d (1+ (sqrt (1+ n))))) 987 ;; (complex (log1+ (/ n d)) (atan i (1+ r)))) 988 ;; (log (1+ x))))) 989 ;; (t 990 (if (and (typep x 'ratio) 991 (< 0.5 x 2)) 992 (setq x (%shortfloat x))) 993 (let ((y (1+ x))) 994 (if (eql y x) 995 (loge y) 996 (let ((e (1 y))) 997 (if (zerop e) 998 (* x 1.0) 999 ( (loge y) (/ ( e x) y))))))) 1000 1001 ;;; helper for complex log 1002 ;;; uses more careful approach when (abs x) is near 1 1003 (defun logabs (x) 1004 (let ((a (abs x))) 1005 (if (< 0.5 a 3) 1006 (let* ((r (realpart x)) 1007 (i (imagpart x)) 1008 (n (if (> (abs r) (abs i)) 1009 (+ (* (1+ r) (1 r)) (* i i)) 1010 (+ (* r r) (* (1+ i) (1 i)))))) 1011 (log1+ (/ n (1+ a)))) 1012 (loge a)))) 1013 1014 (defun %rationalcomplexlog (x prototype &aux ra ia) 1015 (let* ((rx (realpart x)) 1016 (ix (imagpart x)) 1017 (x (abs rx)) 1018 (y (abs ix))) 1019 (if (> y x) 1020 (let ((r (float (/ rx y) 1.0d0))) 1021 (setq ra (+ (%rationallog y 1.0d0) 1022 (/ (log1+ (* r r)) 2))) 1023 (setq ia (atan (if (minusp ix) 1.0d0 1.0d0) r))) 1024 (let ((r (float (/ ix x) 1.0d0))) 1025 (setq ra (+ (%rationallog x 1.0d0) 1026 (/ (log1+ (* r r)) 2))) 1027 (setq ia (atan r (if (minusp rx) 1.0d0 1.0d0))))) 1028 (if (typep prototype 'shortfloat) 1029 (complex (%shortfloat ra) (%shortfloat ia)) 1030 (complex ra ia)))) 807 1031 808 1032 (defun exp (x) … … 854 1078 (cond ((zerop x) x) 855 1079 ((complexp x) 856 (let* ((i (imagpart x))) 857 (if (zerop i) ; has to be a float 858 (let* ((zero (if (typep i 'doublefloat) 0.0d0 0.0f0)) 859 (r (realpart x))) 860 (if (< r zero) 861 (%makecomplex zero (floatsign i (sqrt ( r)))) 862 (%makecomplex (abs (sqrt r)) (floatsign i zero)))) 863 (* (sqrt (abs x)) (cis (/ (phase x) 2)))))) 1080 (let ((rx (realpart x)) 1081 (ix (imagpart x))) 1082 (cond ((rationalp rx) 1083 (if (zerop rx) 1084 (let ((s (sqrt (/ (abs ix) 2)))) 1085 (complex s (if (minusp ix) ( s) s))) 1086 (let* ((s (+ (* rx rx) (* ix ix))) 1087 (d (if (ratiop s) 1088 (/ (isqrt (%numerator s)) 1089 (isqrt (%denominator s))) 1090 (isqrt s)))) 1091 (unless (eql s (* d d)) 1092 (setf d (%doublefloathypot (%doublefloat rx) 1093 (%doublefloat ix)))) 1094 (cond ((minusp rx) 1095 (setq b (sqrt (/ ( d rx) 2))) 1096 (when (minusp ix) 1097 (setq b ( b))) 1098 (setq a (/ ix (+ b b)))) 1099 (t 1100 (setq a (sqrt (/ (+ rx d) 2))) 1101 (setq b (/ ix (+ a a))))) 1102 (if (rationalp a) 1103 (complex a b) 1104 (complex (%shortfloat a) (%shortfloat b)))))) 1105 ((minusp rx) 1106 (if (zerop ix) 1107 (complex 0 (floatsign ix (sqrt ( rx)))) 1108 (let ((shift (cond ((< rx 1) 3) 1109 ((and (> rx 5.9604645E8) (< (abs ix) 5.9604645E8)) 25) 1110 (t 1)))) 1111 (setq rx (scalefloat rx shift)) 1112 (let ((s (fsqrt ( (abs (complex rx (scalefloat ix shift))) rx)))) 1113 (setq b (scalefloat s (ash ( 1 shift) 1))) 1114 (when (minusp ix) 1115 (setq b ( b))) 1116 (setq a (/ ix (scalefloat b 1))) 1117 (complex a b))))) 1118 (t 1119 (if (zerop ix) 1120 (complex (sqrt rx) ix) 1121 (let ((shift (cond ((> rx 1) 3) 1122 ((and (< rx 5.9604645E8) (< (abs ix) 5.9604645E8)) 25) 1123 (t 1)))) 1124 (setq rx (scalefloat rx shift)) 1125 (let ((s (fsqrt (+ rx (abs (complex rx (scalefloat ix shift))))))) 1126 (setq a (scalefloat s (ash ( 1 shift) 1))) 1127 (setq b (/ ix (scalefloat a 1))) 1128 (complex a b)))))))) 864 1129 ((minusp x) (complex 0 (sqrt ( x)))) 865 1130 ((floatp x) … … 874 1139 (/ a b)) 875 1140 (t 876 #+32bittarget 877 (target::withstackshortfloats ((f1)) 878 (fsqrt (%shortfloat x f1))) 879 #+64bittarget 880 (fsqrt (%shortfloat x))))) 1141 (float (fsqrt (float x 0.0d0)) 1.0s0)))) 881 1142 882 1143
Note: See TracChangeset
for help on using the changeset viewer.