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

Propagate r15683 to trunk.

File:
1 edited

Legend:

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

    r15640 r15685  
    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
Note: See TracChangeset for help on using the changeset viewer.