Changeset 15683


Ignore:
Timestamp:
Feb 4, 2013, 6:39:52 PM (6 years ago)
Author:
gb
Message:

Implement the rest of Dave Findlay's "reference" math library.
Fixes ticket:869 in the trunk.

Location:
trunk/source
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/level-0/l0-float.lisp

    r15646 r15683  
    2323  (require "NUMBER-MACROS")
    2424  (require :number-case-macro)
     25  (defconstant two^23 (ash 1 23))
    2526  (defconstant single-float-pi (coerce pi 'single-float))
    2627  (defconstant double-float-half-pi (/ pi 2))
     
    693694(defun sin (x)
    694695  "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 'double-float)
    701       (%double-float-sin! x (%make-dfloat))
    702       #+32-bit-target
    703       (target::with-stack-short-floats ((sx x))
    704         (%single-float-sin! sx (%make-sfloat)))
    705       #+64-bit-target
    706       (%single-float-sin (%short-float 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 'double-float)
     704           (imagpart (%extended-cis x))
     705           (%short-float (imagpart (%extended-cis x)))))
     706        ((typep x 'double-float)
     707         (%double-float-sin! x (%make-dfloat)))
     708        (t
     709         #+32-bit-target
     710         (target::with-stack-short-floats ((sx x))
     711           (%single-float-sin! sx (%make-sfloat)))
     712         #+64-bit-target
     713         (%single-float-sin (%short-float x)))))
    707714
    708715(defun cos (x)
    709716  "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 'double-float)
    716       (%double-float-cos! x (%make-dfloat))
    717       #+32-bit-target
    718       (target::with-stack-short-floats ((sx x))
    719         (%single-float-cos! sx (%make-sfloat)))
    720       #+64-bit-target
    721       (%single-float-cos (%short-float 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 'double-float)
     725           (realpart (%extended-cis x))
     726           (%short-float (realpart (%extended-cis x)))))
     727        ((typep x 'double-float)
     728         (%double-float-cos! x (%make-dfloat)))
     729        (t
     730         #+32-bit-target
     731         (target::with-stack-short-floats ((sx x))
     732           (%single-float-cos! sx (%make-sfloat)))
     733         #+64-bit-target
     734         (%single-float-cos (%short-float x)))))
    722735
    723736(defun tan (x)
    724737  "Return the tangent of NUMBER."
    725   (if (complexp x)
    726     (/ (sin x) (cos x))
    727     (if (typep x 'double-float)
    728       (%double-float-tan! x (%make-dfloat))
    729       #+32-bit-target
    730       (target::with-stack-short-floats ((sx x))
    731         (%single-float-tan! sx (%make-sfloat)))
    732       #+64-bit-target
    733       (%single-float-tan (%short-float 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 (%extended-cis x)))
     756           (if (typep x 'double-float)
     757             (/ (imagpart c) (realpart c))
     758             (%short-float (/ (imagpart c) (realpart c))))))
     759        ((typep x 'double-float)
     760         (%double-float-tan! x (%make-dfloat)))
     761        (t
     762         #+32-bit-target
     763         (target::with-stack-short-floats ((sx x))
     764           (%single-float-tan! sx (%make-sfloat)))
     765         #+64-bit-target
     766         (%single-float-tan (%short-float 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 pi-vector
     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 %extended-cis (x)
     797  (let (size pi-size)
     798    (typecase x
     799      (integer (setq size (1- (integer-length (abs x)))))
     800      (ratio (setq size (- (integer-length (abs (numerator x)))
     801                           (integer-length (denominator x)))))
     802      (float (multiple-value-bind (mantissa exponent sign)
     803                                  (integer-decode-float x)
     804               (setq size (+ (1- (integer-length mantissa)) exponent))
     805               (setq x (* sign mantissa (expt 2 exponent))))))
     806    (setq pi-size (ceiling (+ size 64) 100))
     807    (cond ((< pi-size 1) (setq pi-size 1))
     808          ((> pi-size 17) (setq pi-size 17)))
     809    (let* ((2pi-approx (/ (aref pi-vector (1- pi-size))
     810                          (ash 1 (1- (* 100 pi-size)))))
     811           (reduced-x (rem x 2pi-approx))
     812           (x0 (float reduced-x 1.0d0))
     813           (x1 (- reduced-x (rational x0))))
     814      (* (cis x0) (cis (float x1 1.0d0))))))
    735815
    736816
     
    10461126  (> (realpart n) 0))
    10471127
     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;;;
    10481135(defun expt (b e)
    10491136  "Return BASE raised to the POWER."
     
    10521139         (if (minusp e) (/ 1 (%integer-power b (- e))) (%integer-power b e)))
    10531140        ((zerop b)
    1054          (if (plusp (realpart e)) b (report-bad-arg e '(satisfies positive-realpart-p))))
    1055         ((and (realp b) (plusp b) (realp e))
     1141         (if (plusp (realpart e)) (* b e) (report-bad-arg 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 most-positive-short-float))
     1147                (ratio (cond ((< b 0.5)
     1148                              (>= b least-positive-normalized-short-float))
     1149                             ((> b 3)
     1150                              (<= b most-positive-short-float))))
     1151                (t t)))
     1152         ;; assumes that the library routines are accurate
     1153         ;; (if not, just excise this whole clause)
    10561154         (if (or (typep b 'double-float)
    10571155                 (typep e 'double-float))
     
    10611159           #+32-bit-target
    10621160           (target::with-stack-short-floats ((b1 b)
    1063                                      (e1 e))
     1161                                             (e1 e))
    10641162             (%single-float-expt! b1 e1 (%make-sfloat)))
    10651163           #+64-bit-target
    10661164           (%single-float-expt (%short-float b) (%short-float e))
    10671165           ))
    1068         ((typep (realpart e) 'double-float)
    1069          ;; Avoid intermediate single-float result from LOG
    1070          (let ((promoted-base (* 1d0 b)))
    1071            (exp (* e (log promoted-base)))))
    1072         (t (exp (* e (log b))))))
     1166        ((typep b 'rational)
     1167         (let ((r (exp (* e (%rational-log b 1.0d0)))))
     1168           (cond ((typep (realpart e) 'double-float)
     1169                  r)
     1170                 ((typep r 'complex)
     1171                  (complex (%short-float (realpart r)) (%short-float (imagpart r))))
     1172                 (t
     1173                  (%short-float r)))))
     1174        ((typep (realpart b) 'rational)
     1175         (let ((r (exp (* e (%rational-complex-log b 1.0d0)))))
     1176           (if (typep (realpart e) 'double-float)
     1177             r
     1178             (complex (%short-float (realpart r)) (%short-float (imagpart r))))))
     1179        (t
     1180         ;; type upgrade b without losing -0.0 ...
     1181         (let ((r (exp (* e (log-e (* b 1.0d0))))))
     1182           (cond ((or (typep (realpart b) 'double-float)
     1183                      (typep (realpart e) 'double-float))
     1184                  r)
     1185                 ((typep r 'complex)
     1186                  (complex (%short-float (realpart r)) (%short-float (imagpart r))))
     1187                 (t
     1188                  (%short-float r)))))))
    10731189
    10741190
     
    11451261(defun asin (x)
    11461262  "Return the arc sine of NUMBER."
    1147   (number-case x
    1148     (complex
    1149       (let ((sqrt-1-x (sqrt (- 1 x)))
    1150             (sqrt-1+x (sqrt (+ 1 x))))
    1151         (complex (atan (realpart x)
    1152                        (realpart (* sqrt-1-x sqrt-1+x)))
    1153                  (asinh (imagpart (* (conjugate sqrt-1-x)
    1154                                      sqrt-1+x))))))
    1155     (double-float
    1156      (locally (declare (type double-float x))
    1157        (if (and (<= -1.0d0 x)
    1158                 (<= x 1.0d0))
    1159          (%double-float-asin! x (%make-dfloat))
    1160          (let* ((temp (+ (complex -0.0d0 x)
    1161                          (sqrt (- 1.0d0 (the double-float (* x x)))))))
    1162            (complex (phase temp) (- (log (abs temp))))))))
    1163     ((short-float rational)
    1164      #+32-bit-target
    1165      (let* ((x1 (%make-sfloat)))
    1166        (declare (dynamic-extent x1))
    1167        (if (and (realp x)
    1168                 (<= -1.0s0 (setq x (%short-float x x1)))
    1169                 (<= x 1.0s0))
    1170          (%single-float-asin! x1 (%make-sfloat))
    1171          (progn
    1172            (setq x (+ (complex (- (imagpart x)) (realpart x))
    1173                       (sqrt (- 1 (* x x)))))
    1174            (complex (phase x) (- (log (abs x)))))))
    1175      #+64-bit-target
    1176      (if (and (realp x)
    1177               (<= -1.0s0 (setq x (%short-float x)))
    1178               (<= x 1.0s0))
    1179          (%single-float-asin 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 'double-float)
     1264              (locally (declare (type double-float x))
     1265                (and (<= -1.0d0 x)
     1266                     (<= x 1.0d0))))
     1267         (%double-float-asin! x (%make-dfloat)))
     1268        ((and (typep x 'single-float)
     1269              (locally (declare (type single-float x))
     1270                (and (<= -1.0s0 x)
     1271                     (<= x 1.0s0))))
     1272         #+32-bit-target
     1273         (%single-float-asin! x (%make-sfloat))
     1274         #+64-bit-target
     1275         (%single-float-asin 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 single-float-half-pi)
     1282             (-1 #.(- single-float-half-pi)))
     1283           (atan x (sqrt (- 1 (* x x))))))
     1284        (t
     1285         (%complex-asin/acos x nil))
     1286        ))
    11851287
    11861288
    11871289(defun acos (x)
    11881290  "Return the arc cosine of NUMBER."
    1189   (number-case x
    1190     (complex
    1191      (if (typep (realpart x) 'double-float)
    1192        (- double-float-half-pi (asin x))
    1193        (- single-float-half-pi (asin x))))
    1194     (double-float
    1195      (locally (declare (type double-float x))
    1196        (if (and (<= -1.0d0 x)
    1197                 (<= x 1.0d0))
    1198          (%double-float-acos! x (%make-dfloat))
    1199          (- double-float-half-pi (asin x)))))
    1200     ((short-float rational)
    1201      #+32-bit-target
    1202      (target::with-stack-short-floats ((sx x))
    1203         (locally
    1204             (declare (type short-float sx))
    1205           (if (and (<= -1.0s0 sx)
    1206                    (<= sx 1.0s0))
    1207             (%single-float-acos! sx (%make-sfloat))
    1208             (- single-float-half-pi (asin sx)))))
    1209      #+64-bit-target
    1210      (let* ((sx (%short-float x)))
    1211        (declare (type short-float sx))
    1212        (if (and (<= -1.0s0 sx)
    1213                 (<= sx 1.0s0))
    1214          (%single-float-acos sx)
    1215          (- single-float-half-pi (asin sx))))
    1216      )))
     1291  (cond ((and (typep x 'double-float)
     1292              (locally (declare (type double-float x))
     1293                (and (<= -1.0d0 x)
     1294                     (<= x 1.0d0))))
     1295         (%double-float-acos! x (%make-dfloat)))
     1296        ((and (typep x 'short-float)
     1297              (locally (declare (type short-float x))
     1298                (and (<= -1.0s0 x)
     1299                     (<= x 1.0s0))))
     1300         #+32-bit-target
     1301         (%single-float-acos! x (%make-sfloat))
     1302         #+64-bit-target
     1303         (%single-float-acos x))
     1304        ((and (typep x 'rational)
     1305              (<= (abs x) 1))
     1306         (if (integerp x)
     1307           (case x
     1308             (0 single-float-half-pi)
     1309             (1 0.0s0)                          ; or simply 0 ??
     1310             (-1 single-float-pi))
     1311           (atan (sqrt (- 1 (* x x))) x)))
     1312        (t
     1313         (%complex-asin/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, 299-335 (Sept. 1997)
     1320(defun %complex-asin/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 ((log-s 0))
     1329        (if (typep m 'double-float)
     1330          (if (> m #.(/ most-positive-double-float 2))
     1331            (setq log-s double-float-log2)
     1332            (setq z (* 2 z)))
     1333          (if (> m #.(/ most-positive-short-float 2))
     1334            (setq log-s single-float-log2)
     1335            (setq z (* 2 z))))
     1336        (if acos
     1337          (i* (+ log-s (log-e z))
     1338              (if (minusp ix) +1 -1))
     1339          (if (minusp ix)
     1340            (i* (+ log-s (log-e (i* z 1))) -1)
     1341            (i* (+ log-s (log-e (i* z -1))) 1))))
     1342      (let ((qrx rx)
     1343            (qx x)
     1344            x-1 y2 s)
     1345        (cond ((rationalp rx)
     1346               (setq x-1 (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 x-1)
     1352                              y)
     1353                             ((> y x-1)
     1354                              (let ((c (/ x-1 y)))
     1355                                (* y (sqrt (1+ (* c c))))))
     1356                             (t
     1357                              (let ((c (/ y x-1)))
     1358                                (* x-1 (sqrt (1+ (* c c)))))))))
     1359              (t
     1360               (setq x-1 (abs (- x 1)))
     1361               (setq y2 (* y y))
     1362               (setq s (if (zerop x-1)
     1363                         y
     1364                         (sqrt (+ (* x-1 x-1) 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+x-1 (+ s x-1))
     1373                            (a+x (+ a x))
     1374                            (ry (if (<= qx 1)
     1375                                  (let ((aa (+ (/ y2 r+x+1) s+x-1)))
     1376                                    (sqrt (/ (* a+x aa) 2)))
     1377                                  (let ((aa (+ (/ a+x r+x+1) (/ a+x s+x-1))))
     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+x-1 (+ s x-1))
     1383                            (ll (if (< qx 1)
     1384                                  (let* ((aa (/ (+ (/ 1 r+x+1) (/ 1 s+x-1)) 2)))
     1385                                    (+ (* aa y2) (* y (sqrt (* aa (1+ a))))))
     1386                                  (let* ((a-1 (/ (+ (/ y2 r+x+1) s+x-1) 2)))
     1387                                    (+ a-1 (sqrt (* a-1 (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 (float-sign (- 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 (float-sign ix ia))
     1402                (if (minusp ix)
     1403                  (setq ia (- ia))))
     1404              (if (> qrx 1)
     1405                (setq ia (- ia)))))
     1406          (complex ra ia))))))
    12171407
    12181408
  • trunk/source/level-0/l0-numbers.lisp

    r15649 r15683  
    15911591(defun cis (theta)
    15921592  "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 'double-float)
     1598           (%extended-cis theta)
     1599           (coerce (%extended-cis theta) '(complex single-float))))
     1600        (t
     1601         (complex (cos theta) (sin theta)))))
    15961602
    15971603
  • trunk/source/level-1/l1-numbers.lisp

    r15648 r15683  
    734734    temp))
    735735
     736;;; On current (early 2013) versions of x8664 Linux,
     737;;; (#_atan2 most-positive-double-float most-positive-double-float)
     738;;; returns the correct answer but generates an intermediate
     739;;; invalid-operation exception. #_atan2[f] is documented to never
     740;;; raise fp exceptions, so don't check for them.
    736741(defun %double-float-atan2! (x y result)
    737742  (declare (double-float x y result))
    738   (with-stack-double-floats ((temp))
    739     #+arm-target (%set-fpscr-status 0)
    740     (%setf-double-float TEMP (#_atan2 x y))
    741     (%df-check-exception-2 'atan2 x y (%ffi-exception-status))
    742     (%setf-double-float result TEMP)))
     743  (%setf-double-float result (#_atan2 x y)))
    743744
    744745#+32-bit-target
    745746(defun %single-float-atan2! (x y result)
    746747  (declare (single-float x y result))
    747   (target::with-stack-short-floats ((temp))
    748     #+arm-target (%set-fpscr-status 0)
    749     (%setf-short-float TEMP (#_atan2f x y))
    750     (%sf-check-exception-2 'atan2 x y (%ffi-exception-status))
    751     (%setf-short-float result TEMP)))
     748  (%setf-short-float result (#_atan2f x y)))
    752749
    753750#+64-bit-target
    754751(defun %single-float-atan2 (x y)
    755752  (declare (single-float x y))
    756   (let* ((result (#_atan2f x y)))
    757     (%sf-check-exception-2 'atan2 x y (%ffi-exception-status))
    758     result))
     753  (#_atan2f x y))
    759754
    760755(defun %double-float-exp! (n result)
  • trunk/source/lib/numbers.lisp

    r14665 r15683  
    641641(defun signum (x)
    642642  "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)))))
    644654        ((rationalp x) (if (plusp x) 1 (if (zerop x) 0 -1)))
    645655        ((zerop x) (float 0.0 x))
     
    678688  "Return the hyperbolic sine of NUMBER."
    679689  (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))))
    681694    (if (typep x 'double-float)
    682695      (%double-float-sinh! x (%make-dfloat))
     
    691704  "Return the hyperbolic cosine of NUMBER."
    692705  (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))))
    694710    (if (typep x 'double-float)
    695711      (%double-float-cosh! x (%make-dfloat))
     
    702718(defun tanh (x)
    703719  "Return the hyperbolic tangent of NUMBER."
    704   (if (complexp x)
    705     (/ (sinh x) (cosh x))
    706     (if (typep x 'double-float)
    707       (%double-float-tanh! x (%make-dfloat))
    708       #+32-bit-target
    709       (target::with-stack-short-floats ((sx x))
    710         (%single-float-tanh! sx (%make-sfloat)))
    711       #+64-bit-target
    712       (%single-float-tanh (%short-float 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 'double-float)
     736         (%double-float-tanh! x (%make-dfloat)))
     737        ((and (typep x 'rational)
     738              (> (abs x) 12))
     739         (if (plusp x) 1.0s0 -1.0s0))
     740        (t
     741         #+32-bit-target
     742         (target::with-stack-short-floats ((sx x))
     743           (%single-float-tanh! sx (%make-sfloat)))
     744         #+64-bit-target
     745         (%single-float-tanh (%short-float x)))))
    713746
    714747(defun asinh (x)
    715748  "Return the hyperbolic arc sine of NUMBER."
    716   (if (complexp x)
    717     (log (+ x (sqrt (+ 1 (* x x)))))
    718     (if (typep x 'double-float)
    719       (%double-float-asinh! x (%make-dfloat))
    720       #+32-bit-target
    721       (target::with-stack-short-floats ((sx x))
    722         (%single-float-asinh! sx (%make-sfloat)))
    723       #+64-bit-target
    724       (%single-float-asinh (%short-float x)))))
    725 
     749  (cond ((typep x 'double-float)
     750         (%double-float-asinh! x (%make-dfloat)))
     751        ((typep x 'short-float)
     752         #+32-bit-target
     753         (%single-float-asinh! x (%make-sfloat))
     754         #+64-bit-target
     755         (%single-float-asinh (%short-float x)))
     756        ((typep x 'rational)
     757         (if (< (abs x) most-positive-short-float)
     758           #+32-bit-target
     759           (target::with-stack-short-floats ((sx x))
     760             (%single-float-asinh! sx (%make-sfloat)))
     761           #+64-bit-target
     762           (%single-float-asinh (%short-float x))
     763           (* (signum x) (log-e (* 2 (abs x))))))
     764        (t
     765         (i* (%complex-asin/acos (i* x) nil) -1))))
     766
     767;;; for complex case, use acos and post-fix the branch cut
    726768(defun acosh (x)
    727769  "Return the hyperbolic arc cosine of NUMBER."
    728   (if (and (realp x) (<= 1.0 x))
    729     (if (typep x 'double-float)
    730       (%double-float-acosh! x (%make-dfloat))
    731       #+32-bit-target
    732       (target::with-stack-short-floats ((sx x))
    733         (%single-float-acosh! sx (%make-sfloat)))
    734       #+64-bit-target
    735       (%single-float-acosh (%short-float x)))
    736     (* 2 (log (+ (sqrt (/ (1+ x) 2)) (sqrt (/ (1- x) 2)))))))
     770  (cond ((and (typep x 'double-float)
     771              (locally (declare (type double-float x))
     772                (<= 1.0d0 x)))
     773         (%double-float-acosh! x (%make-dfloat)))
     774        ((and (typep x 'short-float)
     775              (locally (declare (type short-float x))
     776                (<= 1.0s0 x)))
     777         #+32-bit-target
     778         (%single-float-acosh! x (%make-sfloat))
     779         #+64-bit-target
     780         (%single-float-acosh (%short-float x)))
     781        ((and (typep x 'rational)
     782              (<= 1 x))
     783         (cond ((< x 2)
     784                (log1+ (+ (- x 1) (sqrt (- (* x x) 1)))))
     785               ((<= x most-positive-short-float)
     786                #+32-bit-target
     787                (target::with-stack-short-floats ((x1 x))
     788                  (%single-float-acosh! x1 (%make-sfloat)))
     789                #+64-bit-target
     790                (%single-float-acosh (%short-float x)))
     791               (t
     792                (log-e (* 2 x)))))
     793        (t
     794         (let ((sign (and (typep x 'complex)
     795                          (let ((ix (imagpart x)))
     796                            (typecase ix
     797                              (double-float (%double-float-sign ix))
     798                              (single-float (%short-float-sign ix))
     799                              (t (minusp ix)))))))
     800           (i* (%complex-asin/acos x t) (if sign -1 1))))))
    737801
    738802(defun atanh (x)
    739803  "Return the hyperbolic arc tangent of NUMBER."
    740   (if (and (realp x) (<= -1.0 (setq x (float x)) 1.0))
    741     (if (typep x 'double-float)
    742       (%double-float-atanh! x (%make-dfloat))
    743       #+32-bit-target
    744       (%single-float-atanh! x (%make-sfloat))
    745       #+64-bit-target
    746       (%single-float-atanh x))
    747     (/ (- (log (+ 1 x)) (log (- 1 x))) 2)))
     804  (cond ((and (typep x 'double-float)
     805              (locally (declare (type double-float x))
     806                (and (<= -1.0d0 x)
     807                     (<= x 1.0d0))))
     808         (%double-float-atanh! x (%make-dfloat)))
     809        ((and (typep x 'short-float)
     810              (locally (declare (type short-float x))
     811                (and (<= -1.0s0 x)
     812                     (<= x 1.0s0))))
     813         #+32-bit-target
     814         (%single-float-atanh! x (%make-sfloat))
     815         #+64-bit-target
     816         (%single-float-atanh x))
     817        ((and (typep x 'rational)
     818              (<= (abs x) 1))
     819         (let ((n (numerator x))
     820               (d (denominator x)))
     821           (/ (log-e (/ (+ d n) (- d n))) 2)))
     822        (t
     823         (let ((r (realpart x)))
     824           (if (zerop r)
     825             (complex r (atan (imagpart x)))
     826             (%complex-atanh x))))))
    748827
    749828(defun ffloor (number &optional divisor)
Note: See TracChangeset for help on using the changeset viewer.