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/lib/numbers.lisp

    r14665 r15685  
    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.