Changeset 15683
 Timestamp:
 Feb 4, 2013, 6:39:52 PM (6 years ago)
 Location:
 trunk/source
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

trunk/source/level0/l0float.lisp
r15646 r15683 23 23 (require "NUMBERMACROS") 24 24 (require :numbercasemacro) 25 (defconstant two^23 (ash 1 23)) 25 26 (defconstant singlefloatpi (coerce pi 'singlefloat)) 26 27 (defconstant doublefloathalfpi (/ pi 2)) … … 693 694 (defun sin (x) 694 695 "Return the sine of NUMBER." 695 (if (complexp x) 696 (let* ((r (realpart x)) 697 (i (imagpart x))) 698 (complex (* (sin r) (cosh i)) 699 (* (cos r) (sinh i)))) 700 (if (typep x 'doublefloat) 701 (%doublefloatsin! x (%makedfloat)) 702 #+32bittarget 703 (target::withstackshortfloats ((sx x)) 704 (%singlefloatsin! sx (%makesfloat))) 705 #+64bittarget 706 (%singlefloatsin (%shortfloat x))))) 696 (cond ((complexp x) 697 (let* ((r (realpart x)) 698 (i (imagpart x))) 699 (complex (* (sin r) (cosh i)) 700 (* (cos r) (sinh i))))) 701 ((or (typep x 'ratio) 702 (> (abs x) two^23)) 703 (if (typep x 'doublefloat) 704 (imagpart (%extendedcis x)) 705 (%shortfloat (imagpart (%extendedcis x))))) 706 ((typep x 'doublefloat) 707 (%doublefloatsin! x (%makedfloat))) 708 (t 709 #+32bittarget 710 (target::withstackshortfloats ((sx x)) 711 (%singlefloatsin! sx (%makesfloat))) 712 #+64bittarget 713 (%singlefloatsin (%shortfloat x))))) 707 714 708 715 (defun cos (x) 709 716 "Return the cosine of NUMBER." 710 (if (complexp x) 711 (let* ((r (realpart x)) 712 (i (imagpart x))) 713 (complex (* (cos r) (cosh i)) 714 ( (* (sin r) (sinh i))))) 715 (if (typep x 'doublefloat) 716 (%doublefloatcos! x (%makedfloat)) 717 #+32bittarget 718 (target::withstackshortfloats ((sx x)) 719 (%singlefloatcos! sx (%makesfloat))) 720 #+64bittarget 721 (%singlefloatcos (%shortfloat x))))) 717 (cond ((complexp x) 718 (let* ((r (realpart x)) 719 (i (imagpart x))) 720 (complex (* (cos r) (cosh i)) 721 ( (* (sin r) (sinh i)))))) 722 ((or (typep x 'ratio) 723 (> (abs x) two^23)) 724 (if (typep x 'doublefloat) 725 (realpart (%extendedcis x)) 726 (%shortfloat (realpart (%extendedcis x))))) 727 ((typep x 'doublefloat) 728 (%doublefloatcos! x (%makedfloat))) 729 (t 730 #+32bittarget 731 (target::withstackshortfloats ((sx x)) 732 (%singlefloatcos! sx (%makesfloat))) 733 #+64bittarget 734 (%singlefloatcos (%shortfloat x))))) 722 735 723 736 (defun tan (x) 724 737 "Return the tangent of NUMBER." 725 (if (complexp x) 726 (/ (sin x) (cos x)) 727 (if (typep x 'doublefloat) 728 (%doublefloattan! x (%makedfloat)) 729 #+32bittarget 730 (target::withstackshortfloats ((sx x)) 731 (%singlefloattan! sx (%makesfloat))) 732 #+64bittarget 733 (%singlefloattan (%shortfloat x)) 734 ))) 738 (cond ((complexp x) 739 (let ((r (realpart x)) 740 (i (imagpart x))) 741 (if (zerop i) 742 (complex (tan r) i) 743 (let* ((tx (tan r)) 744 (ty (tanh i)) 745 (tx2 (* tx tx)) 746 (d (1+ (* tx2 (* ty ty)))) 747 (n (if (> (abs i) 20) 748 (* 4 (exp (* 2 (abs i)))) 749 (let ((c (cosh i))) 750 (/ (* c c)))))) 751 (complex (/ (* n tx) d) 752 (/ (* ty (1+ tx2)) d)))))) 753 ((or (typep x 'ratio) 754 (> (abs x) two^23)) 755 (let ((c (%extendedcis x))) 756 (if (typep x 'doublefloat) 757 (/ (imagpart c) (realpart c)) 758 (%shortfloat (/ (imagpart c) (realpart c)))))) 759 ((typep x 'doublefloat) 760 (%doublefloattan! x (%makedfloat))) 761 (t 762 #+32bittarget 763 (target::withstackshortfloats ((sx x)) 764 (%singlefloattan! sx (%makesfloat))) 765 #+64bittarget 766 (%singlefloattan (%shortfloat x)) 767 ))) 768 769 770 ;;; Helper function for sin/cos/tan for ratio or large arguments 771 ;;; (Will become inaccurate for ridiculously large arguments though) 772 ;;; Does not assume that float library returns accurate values for large arguments 773 ;;; (many don't) 774 775 ;;; hexdecimal representations of pi at various precisions 776 (defconstant pivector 777 #(#x3243F6A8885A308D313198A2E0 778 #x3243F6A8885A308D313198A2E03707344A4093822299F31D008 779 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D 780 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC 781 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B5470 782 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310B 783 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB 784 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045 785 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70 786 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871 787 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D9 788 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD588215 789 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D 790 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085 791 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB 792 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3E 793 #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3ED71577C1BD314B2778AF2FDA5 794 )) 795 796 (defun %extendedcis (x) 797 (let (size pisize) 798 (typecase x 799 (integer (setq size (1 (integerlength (abs x))))) 800 (ratio (setq size ( (integerlength (abs (numerator x))) 801 (integerlength (denominator x))))) 802 (float (multiplevaluebind (mantissa exponent sign) 803 (integerdecodefloat x) 804 (setq size (+ (1 (integerlength mantissa)) exponent)) 805 (setq x (* sign mantissa (expt 2 exponent)))))) 806 (setq pisize (ceiling (+ size 64) 100)) 807 (cond ((< pisize 1) (setq pisize 1)) 808 ((> pisize 17) (setq pisize 17))) 809 (let* ((2piapprox (/ (aref pivector (1 pisize)) 810 (ash 1 (1 (* 100 pisize))))) 811 (reducedx (rem x 2piapprox)) 812 (x0 (float reducedx 1.0d0)) 813 (x1 ( reducedx (rational x0)))) 814 (* (cis x0) (cis (float x1 1.0d0)))))) 735 815 736 816 … … 1046 1126 (> (realpart n) 0)) 1047 1127 1128 ;;; (It may be possible to do something with rational exponents, e.g. so that 1129 ;;; (expt x 1/2) == (sqrt x) 1130 ;;; (expt x 3/2) == (expt (sqrt x) 3) ** NOT (sqrt (expt x 3)) !! ** 1131 ;;; or (* x (sqrt x)) 1132 ;;; (expt x 1/8) == (sqrt (sqrt (sqrt x))) 1133 ;;; even, in the case of rational x, returning a rational result if possible.) 1134 ;;; 1048 1135 (defun expt (b e) 1049 1136 "Return BASE raised to the POWER." … … 1052 1139 (if (minusp e) (/ 1 (%integerpower b ( e))) (%integerpower b e))) 1053 1140 ((zerop b) 1054 (if (plusp (realpart e)) b (reportbadarg e '(satisfies positiverealpartp)))) 1055 ((and (realp b) (plusp b) (realp e)) 1141 (if (plusp (realpart e)) (* b e) (reportbadarg e '(satisfies plusp)))) 1142 ((and (realp b) (plusp b) (realp e) 1143 ; escape out very small or very large rationals 1144 ;  avoid problems converting to float 1145 (typecase b 1146 (bignum (<= b mostpositiveshortfloat)) 1147 (ratio (cond ((< b 0.5) 1148 (>= b leastpositivenormalizedshortfloat)) 1149 ((> b 3) 1150 (<= b mostpositiveshortfloat)))) 1151 (t t))) 1152 ;; assumes that the library routines are accurate 1153 ;; (if not, just excise this whole clause) 1056 1154 (if (or (typep b 'doublefloat) 1057 1155 (typep e 'doublefloat)) … … 1061 1159 #+32bittarget 1062 1160 (target::withstackshortfloats ((b1 b) 1063 (e1 e))1161 (e1 e)) 1064 1162 (%singlefloatexpt! b1 e1 (%makesfloat))) 1065 1163 #+64bittarget 1066 1164 (%singlefloatexpt (%shortfloat b) (%shortfloat e)) 1067 1165 )) 1068 ((typep (realpart e) 'doublefloat) 1069 ;; Avoid intermediate singlefloat result from LOG 1070 (let ((promotedbase (* 1d0 b))) 1071 (exp (* e (log promotedbase))))) 1072 (t (exp (* e (log b)))))) 1166 ((typep b 'rational) 1167 (let ((r (exp (* e (%rationallog b 1.0d0))))) 1168 (cond ((typep (realpart e) 'doublefloat) 1169 r) 1170 ((typep r 'complex) 1171 (complex (%shortfloat (realpart r)) (%shortfloat (imagpart r)))) 1172 (t 1173 (%shortfloat r))))) 1174 ((typep (realpart b) 'rational) 1175 (let ((r (exp (* e (%rationalcomplexlog b 1.0d0))))) 1176 (if (typep (realpart e) 'doublefloat) 1177 r 1178 (complex (%shortfloat (realpart r)) (%shortfloat (imagpart r)))))) 1179 (t 1180 ;; type upgrade b without losing 0.0 ... 1181 (let ((r (exp (* e (loge (* b 1.0d0)))))) 1182 (cond ((or (typep (realpart b) 'doublefloat) 1183 (typep (realpart e) 'doublefloat)) 1184 r) 1185 ((typep r 'complex) 1186 (complex (%shortfloat (realpart r)) (%shortfloat (imagpart r)))) 1187 (t 1188 (%shortfloat r))))))) 1073 1189 1074 1190 … … 1145 1261 (defun asin (x) 1146 1262 "Return the arc sine of NUMBER." 1147 (numbercase x 1148 (complex 1149 (let ((sqrt1x (sqrt ( 1 x))) 1150 (sqrt1+x (sqrt (+ 1 x)))) 1151 (complex (atan (realpart x) 1152 (realpart (* sqrt1x sqrt1+x))) 1153 (asinh (imagpart (* (conjugate sqrt1x) 1154 sqrt1+x)))))) 1155 (doublefloat 1156 (locally (declare (type doublefloat x)) 1157 (if (and (<= 1.0d0 x) 1158 (<= x 1.0d0)) 1159 (%doublefloatasin! x (%makedfloat)) 1160 (let* ((temp (+ (complex 0.0d0 x) 1161 (sqrt ( 1.0d0 (the doublefloat (* x x))))))) 1162 (complex (phase temp) ( (log (abs temp)))))))) 1163 ((shortfloat rational) 1164 #+32bittarget 1165 (let* ((x1 (%makesfloat))) 1166 (declare (dynamicextent x1)) 1167 (if (and (realp x) 1168 (<= 1.0s0 (setq x (%shortfloat x x1))) 1169 (<= x 1.0s0)) 1170 (%singlefloatasin! x1 (%makesfloat)) 1171 (progn 1172 (setq x (+ (complex ( (imagpart x)) (realpart x)) 1173 (sqrt ( 1 (* x x))))) 1174 (complex (phase x) ( (log (abs x))))))) 1175 #+64bittarget 1176 (if (and (realp x) 1177 (<= 1.0s0 (setq x (%shortfloat x))) 1178 (<= x 1.0s0)) 1179 (%singlefloatasin x) 1180 (progn 1181 (setq x (+ (complex ( (imagpart x)) (realpart x)) 1182 (sqrt ( 1 (* x x))))) 1183 (complex (phase x) ( (log (abs x)))))) 1184 ))) 1263 (cond ((and (typep x 'doublefloat) 1264 (locally (declare (type doublefloat x)) 1265 (and (<= 1.0d0 x) 1266 (<= x 1.0d0)))) 1267 (%doublefloatasin! x (%makedfloat))) 1268 ((and (typep x 'singlefloat) 1269 (locally (declare (type singlefloat x)) 1270 (and (<= 1.0s0 x) 1271 (<= x 1.0s0)))) 1272 #+32bittarget 1273 (%singlefloatasin! x (%makesfloat)) 1274 #+64bittarget 1275 (%singlefloatasin x)) 1276 ((and (typep x 'rational) 1277 (<= (abs x) 1)) 1278 (if (integerp x) 1279 (case x 1280 (0 0.0s0) ; or simply 0 ?? 1281 (1 singlefloathalfpi) 1282 (1 #.( singlefloathalfpi))) 1283 (atan x (sqrt ( 1 (* x x)))))) 1284 (t 1285 (%complexasin/acos x nil)) 1286 )) 1185 1287 1186 1288 1187 1289 (defun acos (x) 1188 1290 "Return the arc cosine of NUMBER." 1189 (numbercase x 1190 (complex 1191 (if (typep (realpart x) 'doublefloat) 1192 ( doublefloathalfpi (asin x)) 1193 ( singlefloathalfpi (asin x)))) 1194 (doublefloat 1195 (locally (declare (type doublefloat x)) 1196 (if (and (<= 1.0d0 x) 1197 (<= x 1.0d0)) 1198 (%doublefloatacos! x (%makedfloat)) 1199 ( doublefloathalfpi (asin x))))) 1200 ((shortfloat rational) 1201 #+32bittarget 1202 (target::withstackshortfloats ((sx x)) 1203 (locally 1204 (declare (type shortfloat sx)) 1205 (if (and (<= 1.0s0 sx) 1206 (<= sx 1.0s0)) 1207 (%singlefloatacos! sx (%makesfloat)) 1208 ( singlefloathalfpi (asin sx))))) 1209 #+64bittarget 1210 (let* ((sx (%shortfloat x))) 1211 (declare (type shortfloat sx)) 1212 (if (and (<= 1.0s0 sx) 1213 (<= sx 1.0s0)) 1214 (%singlefloatacos sx) 1215 ( singlefloathalfpi (asin sx)))) 1216 ))) 1291 (cond ((and (typep x 'doublefloat) 1292 (locally (declare (type doublefloat x)) 1293 (and (<= 1.0d0 x) 1294 (<= x 1.0d0)))) 1295 (%doublefloatacos! x (%makedfloat))) 1296 ((and (typep x 'shortfloat) 1297 (locally (declare (type shortfloat x)) 1298 (and (<= 1.0s0 x) 1299 (<= x 1.0s0)))) 1300 #+32bittarget 1301 (%singlefloatacos! x (%makesfloat)) 1302 #+64bittarget 1303 (%singlefloatacos x)) 1304 ((and (typep x 'rational) 1305 (<= (abs x) 1)) 1306 (if (integerp x) 1307 (case x 1308 (0 singlefloathalfpi) 1309 (1 0.0s0) ; or simply 0 ?? 1310 (1 singlefloatpi)) 1311 (atan (sqrt ( 1 (* x x))) x))) 1312 (t 1313 (%complexasin/acos x t)) 1314 )) 1315 1316 ;;; combined complex asin/acos routine 1317 ;;; argument acos is true for acos(z); false for asin(z) 1318 ;;; 1319 ;;; based on Hull, Fairgrieve & Tang, ACM TMS 23, 3, 299335 (Sept. 1997) 1320 (defun %complexasin/acos (z acos) 1321 (let* ((rx (realpart z)) 1322 (ix (imagpart z)) 1323 (x (abs rx)) 1324 (y (abs ix)) 1325 (m (max x y))) 1326 (if (> m 1.8014399E+16) 1327 ;; Large argument escape 1328 (let ((logs 0)) 1329 (if (typep m 'doublefloat) 1330 (if (> m #.(/ mostpositivedoublefloat 2)) 1331 (setq logs doublefloatlog2) 1332 (setq z (* 2 z))) 1333 (if (> m #.(/ mostpositiveshortfloat 2)) 1334 (setq logs singlefloatlog2) 1335 (setq z (* 2 z)))) 1336 (if acos 1337 (i* (+ logs (loge z)) 1338 (if (minusp ix) +1 1)) 1339 (if (minusp ix) 1340 (i* (+ logs (loge (i* z 1))) 1) 1341 (i* (+ logs (loge (i* z 1))) 1)))) 1342 (let ((qrx rx) 1343 (qx x) 1344 x1 y2 s) 1345 (cond ((rationalp rx) 1346 (setq x1 (float (abs ( x 1)))) 1347 (setq rx (float rx)) 1348 (setq x (abs rx)) 1349 (setq y (float y)) 1350 (setq y2 (* y y)) 1351 (setq s (cond ((zerop x1) 1352 y) 1353 ((> y x1) 1354 (let ((c (/ x1 y))) 1355 (* y (sqrt (1+ (* c c)))))) 1356 (t 1357 (let ((c (/ y x1))) 1358 (* x1 (sqrt (1+ (* c c))))))))) 1359 (t 1360 (setq x1 (abs ( x 1))) 1361 (setq y2 (* y y)) 1362 (setq s (if (zerop x1) 1363 y 1364 (sqrt (+ (* x1 x1) y2)))))) 1365 (let* ((x+1 (+ x 1)) 1366 (r (sqrt (+ (* x+1 x+1) y2))) 1367 (a (/ (+ r s) 2)) 1368 (b (/ rx a)) 1369 (ra (if (<= (abs b) 0.6417) 1370 (if acos (acos b) (asin b)) 1371 (let* ((r+x+1 (+ r x+1)) 1372 (s+x1 (+ s x1)) 1373 (a+x (+ a x)) 1374 (ry (if (<= qx 1) 1375 (let ((aa (+ (/ y2 r+x+1) s+x1))) 1376 (sqrt (/ (* a+x aa) 2))) 1377 (let ((aa (+ (/ a+x r+x+1) (/ a+x s+x1)))) 1378 (* y (sqrt (/ aa 2))))))) 1379 (if acos (atan ry rx) (atan rx ry))))) 1380 (ia (if (<= a 1.5) 1381 (let* ((r+x+1 (+ r x+1)) 1382 (s+x1 (+ s x1)) 1383 (ll (if (< qx 1) 1384 (let* ((aa (/ (+ (/ 1 r+x+1) (/ 1 s+x1)) 2))) 1385 (+ (* aa y2) (* y (sqrt (* aa (1+ a)))))) 1386 (let* ((a1 (/ (+ (/ y2 r+x+1) s+x1) 2))) 1387 (+ a1 (sqrt (* a1 (1+ a)))))))) 1388 (log1+ ll)) 1389 (log (+ a (sqrt (1 (* a a)))))))) 1390 ;; final fixup of signs 1391 (if acos 1392 (if (complexp z) 1393 (if (typep ix 'float) 1394 (setq ia (floatsign ( ix) ia)) 1395 (if (plusp ix) 1396 (setq ia ( ia)))) 1397 (if (< qrx 1) 1398 (setq ia ( ia)))) 1399 (if (complexp z) 1400 (if (typep ix 'float) 1401 (setq ia (floatsign ix ia)) 1402 (if (minusp ix) 1403 (setq ia ( ia)))) 1404 (if (> qrx 1) 1405 (setq ia ( ia))))) 1406 (complex ra ia)))))) 1217 1407 1218 1408 
trunk/source/level0/l0numbers.lisp
r15649 r15683 1591 1591 (defun cis (theta) 1592 1592 "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)." 1593 (if (complexp theta) 1594 (error "Argument to CIS is complex: ~S" theta) 1595 (complex (cos theta) (sin theta)))) 1593 (cond ((complexp theta) 1594 (error "Argument to CIS is complex: ~S" theta)) 1595 ((or (typep theta 'ratio) 1596 (> (abs theta) #.(ash 1 23))) 1597 (if (typep theta 'doublefloat) 1598 (%extendedcis theta) 1599 (coerce (%extendedcis theta) '(complex singlefloat)))) 1600 (t 1601 (complex (cos theta) (sin theta))))) 1596 1602 1597 1603 
trunk/source/level1/l1numbers.lisp
r15648 r15683 734 734 temp)) 735 735 736 ;;; On current (early 2013) versions of x8664 Linux, 737 ;;; (#_atan2 mostpositivedoublefloat mostpositivedoublefloat) 738 ;;; returns the correct answer but generates an intermediate 739 ;;; invalidoperation exception. #_atan2[f] is documented to never 740 ;;; raise fp exceptions, so don't check for them. 736 741 (defun %doublefloatatan2! (x y result) 737 742 (declare (doublefloat x y result)) 738 (withstackdoublefloats ((temp)) 739 #+armtarget (%setfpscrstatus 0) 740 (%setfdoublefloat TEMP (#_atan2 x y)) 741 (%dfcheckexception2 'atan2 x y (%ffiexceptionstatus)) 742 (%setfdoublefloat result TEMP))) 743 (%setfdoublefloat result (#_atan2 x y))) 743 744 744 745 #+32bittarget 745 746 (defun %singlefloatatan2! (x y result) 746 747 (declare (singlefloat x y result)) 747 (target::withstackshortfloats ((temp)) 748 #+armtarget (%setfpscrstatus 0) 749 (%setfshortfloat TEMP (#_atan2f x y)) 750 (%sfcheckexception2 'atan2 x y (%ffiexceptionstatus)) 751 (%setfshortfloat result TEMP))) 748 (%setfshortfloat result (#_atan2f x y))) 752 749 753 750 #+64bittarget 754 751 (defun %singlefloatatan2 (x y) 755 752 (declare (singlefloat x y)) 756 (let* ((result (#_atan2f x y))) 757 (%sfcheckexception2 'atan2 x y (%ffiexceptionstatus)) 758 result)) 753 (#_atan2f x y)) 759 754 760 755 (defun %doublefloatexp! (n result) 
trunk/source/lib/numbers.lisp
r14665 r15683 641 641 (defun signum (x) 642 642 "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))." 643 (cond ((complexp x) (if (zerop x) x (/ x (abs x)))) 643 (cond ((complexp x) (if (zerop x) 644 x 645 (let ((m (max (abs (realpart x)) 646 (abs (imagpart x))))) 647 (cond ((rationalp m) 648 ;; rescale to avoid intermediate under/overflow 649 (setq x (/ x m))) 650 ((> m #.(ash 1 23)) 651 ;; ensure no overflow for abs 652 (setq x (/ x 2)))) 653 (/ x (abs x))))) 644 654 ((rationalp x) (if (plusp x) 1 (if (zerop x) 0 1))) 645 655 ((zerop x) (float 0.0 x)) … … 678 688 "Return the hyperbolic sine of NUMBER." 679 689 (if (complexp x) 680 (/ ( (exp x) (exp ( x))) 2) 690 (let ((r (realpart x)) 691 (i (imagpart x))) 692 (complex (* (sinh r) (cos i)) 693 (* (cosh r) (sin i)))) 681 694 (if (typep x 'doublefloat) 682 695 (%doublefloatsinh! x (%makedfloat)) … … 691 704 "Return the hyperbolic cosine of NUMBER." 692 705 (if (complexp x) 693 (/ (+ (exp x) (exp ( x))) 2) 706 (let ((r (realpart x)) 707 (i (imagpart x))) 708 (complex (* (cosh r) (cos i)) 709 (* (sinh r) (sin i)))) 694 710 (if (typep x 'doublefloat) 695 711 (%doublefloatcosh! x (%makedfloat)) … … 702 718 (defun tanh (x) 703 719 "Return the hyperbolic tangent of NUMBER." 704 (if (complexp x) 705 (/ (sinh x) (cosh x)) 706 (if (typep x 'doublefloat) 707 (%doublefloattanh! x (%makedfloat)) 708 #+32bittarget 709 (target::withstackshortfloats ((sx x)) 710 (%singlefloattanh! sx (%makesfloat))) 711 #+64bittarget 712 (%singlefloattanh (%shortfloat x))))) 720 (cond ((complexp x) 721 (let ((r (realpart x)) 722 (i (imagpart x))) 723 (if (zerop r) 724 (complex r (tan i)) 725 (let* ((tx (tanh r)) 726 (ty (tan i)) 727 (ty2 (* ty ty)) 728 (d (1+ (* (* tx tx) ty2))) 729 (n (if (> (abs r) 20) 730 (* 4 (exp ( (* 2 (abs r))))) 731 (let ((c (cosh r))) 732 (/ (* c c)))))) 733 (complex (/ (* tx (1+ ty2)) d) 734 (/ (* n ty) d)))))) 735 ((typep x 'doublefloat) 736 (%doublefloattanh! x (%makedfloat))) 737 ((and (typep x 'rational) 738 (> (abs x) 12)) 739 (if (plusp x) 1.0s0 1.0s0)) 740 (t 741 #+32bittarget 742 (target::withstackshortfloats ((sx x)) 743 (%singlefloattanh! sx (%makesfloat))) 744 #+64bittarget 745 (%singlefloattanh (%shortfloat x))))) 713 746 714 747 (defun asinh (x) 715 748 "Return the hyperbolic arc sine of NUMBER." 716 (if (complexp x) 717 (log (+ x (sqrt (+ 1 (* x x))))) 718 (if (typep x 'doublefloat) 719 (%doublefloatasinh! x (%makedfloat)) 720 #+32bittarget 721 (target::withstackshortfloats ((sx x)) 722 (%singlefloatasinh! sx (%makesfloat))) 723 #+64bittarget 724 (%singlefloatasinh (%shortfloat x))))) 725 749 (cond ((typep x 'doublefloat) 750 (%doublefloatasinh! x (%makedfloat))) 751 ((typep x 'shortfloat) 752 #+32bittarget 753 (%singlefloatasinh! x (%makesfloat)) 754 #+64bittarget 755 (%singlefloatasinh (%shortfloat x))) 756 ((typep x 'rational) 757 (if (< (abs x) mostpositiveshortfloat) 758 #+32bittarget 759 (target::withstackshortfloats ((sx x)) 760 (%singlefloatasinh! sx (%makesfloat))) 761 #+64bittarget 762 (%singlefloatasinh (%shortfloat x)) 763 (* (signum x) (loge (* 2 (abs x)))))) 764 (t 765 (i* (%complexasin/acos (i* x) nil) 1)))) 766 767 ;;; for complex case, use acos and postfix the branch cut 726 768 (defun acosh (x) 727 769 "Return the hyperbolic arc cosine of NUMBER." 728 (if (and (realp x) (<= 1.0 x)) 729 (if (typep x 'doublefloat) 730 (%doublefloatacosh! x (%makedfloat)) 731 #+32bittarget 732 (target::withstackshortfloats ((sx x)) 733 (%singlefloatacosh! sx (%makesfloat))) 734 #+64bittarget 735 (%singlefloatacosh (%shortfloat x))) 736 (* 2 (log (+ (sqrt (/ (1+ x) 2)) (sqrt (/ (1 x) 2))))))) 770 (cond ((and (typep x 'doublefloat) 771 (locally (declare (type doublefloat x)) 772 (<= 1.0d0 x))) 773 (%doublefloatacosh! x (%makedfloat))) 774 ((and (typep x 'shortfloat) 775 (locally (declare (type shortfloat x)) 776 (<= 1.0s0 x))) 777 #+32bittarget 778 (%singlefloatacosh! x (%makesfloat)) 779 #+64bittarget 780 (%singlefloatacosh (%shortfloat x))) 781 ((and (typep x 'rational) 782 (<= 1 x)) 783 (cond ((< x 2) 784 (log1+ (+ ( x 1) (sqrt ( (* x x) 1))))) 785 ((<= x mostpositiveshortfloat) 786 #+32bittarget 787 (target::withstackshortfloats ((x1 x)) 788 (%singlefloatacosh! x1 (%makesfloat))) 789 #+64bittarget 790 (%singlefloatacosh (%shortfloat x))) 791 (t 792 (loge (* 2 x))))) 793 (t 794 (let ((sign (and (typep x 'complex) 795 (let ((ix (imagpart x))) 796 (typecase ix 797 (doublefloat (%doublefloatsign ix)) 798 (singlefloat (%shortfloatsign ix)) 799 (t (minusp ix))))))) 800 (i* (%complexasin/acos x t) (if sign 1 1)))))) 737 801 738 802 (defun atanh (x) 739 803 "Return the hyperbolic arc tangent of NUMBER." 740 (if (and (realp x) (<= 1.0 (setq x (float x)) 1.0)) 741 (if (typep x 'doublefloat) 742 (%doublefloatatanh! x (%makedfloat)) 743 #+32bittarget 744 (%singlefloatatanh! x (%makesfloat)) 745 #+64bittarget 746 (%singlefloatatanh x)) 747 (/ ( (log (+ 1 x)) (log ( 1 x))) 2))) 804 (cond ((and (typep x 'doublefloat) 805 (locally (declare (type doublefloat x)) 806 (and (<= 1.0d0 x) 807 (<= x 1.0d0)))) 808 (%doublefloatatanh! x (%makedfloat))) 809 ((and (typep x 'shortfloat) 810 (locally (declare (type shortfloat x)) 811 (and (<= 1.0s0 x) 812 (<= x 1.0s0)))) 813 #+32bittarget 814 (%singlefloatatanh! x (%makesfloat)) 815 #+64bittarget 816 (%singlefloatatanh x)) 817 ((and (typep x 'rational) 818 (<= (abs x) 1)) 819 (let ((n (numerator x)) 820 (d (denominator x))) 821 (/ (loge (/ (+ d n) ( d n))) 2))) 822 (t 823 (let ((r (realpart x))) 824 (if (zerop r) 825 (complex r (atan (imagpart x))) 826 (%complexatanh x)))))) 748 827 749 828 (defun ffloor (number &optional divisor)
Note: See TracChangeset
for help on using the changeset viewer.