src/code/irrat.lisp

   1 ;;;; This file contains all the irrational functions. (Actually, most
   2 ;;;; of the work is done by calling out to C.)
   3
   4 ;;;; This software is part of the SBCL system. See the README file for
   5 ;;;; more information.
   6 ;;;;
   7 ;;;; This software is derived from the CMU CL system, which was
   8 ;;;; written at Carnegie Mellon University and released into the
   9 ;;;; public domain. The software is in the public domain and is
  10 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
  11 ;;;; files for more information.
  12
  13 (in-package "SB!KERNEL")
  14
  15 (file-comment
  16   "$Header$")
  17 \f
  18 ;;;; miscellaneous constants, utility functions, and macros
  19
  20 (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
  21 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
  22
  23 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
  24 ;;; call, and saves number consing to boot.
  25 ;;;
  26 ;;; FIXME: This should be (EVAL-WHEN (COMPILE-EVAL) (SB!XC:DEFMACRO ..)),
  27 ;;; I think.
  28 (defmacro def-math-rtn (name num-args)
  29   (let ((function (intern (concatenate 'simple-string
  30                                        "%"
  31                                        (string-upcase name)))))
  32     `(progn
  33        (proclaim '(inline ,function))
  34        (let ((sb!int::*rogue-export* "DEF-MATH-RTN"))
  35          (export ',function))
  36        (sb!alien:def-alien-routine (,name ,function) double-float
  37          ,@(let ((results nil))
  38              (dotimes (i num-args (nreverse results))
  39                (push (list (intern (format nil "ARG-~D" i))
  40                            'double-float)
  41                      results)))))))
  42
  43 (eval-when (:compile-toplevel :load-toplevel :execute)
  44
  45 (defun handle-reals (function var)
  46   `((((foreach fixnum single-float bignum ratio))
  47      (coerce (,function (coerce ,var 'double-float)) 'single-float))
  48     ((double-float)
  49      (,function ,var))))
  50
  51 ) ; EVAL-WHEN
  52 \f
  53 ;;;; stubs for the Unix math library
  54
  55 ;;; Please refer to the Unix man pages for details about these routines.
  56
  57 ;;; Trigonometric.
  58 #!-x86 (def-math-rtn "sin" 1)
  59 #!-x86 (def-math-rtn "cos" 1)
  60 #!-x86 (def-math-rtn "tan" 1)
  61 (def-math-rtn "asin" 1)
  62 (def-math-rtn "acos" 1)
  63 #!-x86 (def-math-rtn "atan" 1)
  64 #!-x86 (def-math-rtn "atan2" 2)
  65 (def-math-rtn "sinh" 1)
  66 (def-math-rtn "cosh" 1)
  67 (def-math-rtn "tanh" 1)
  68 (def-math-rtn "asinh" 1)
  69 (def-math-rtn "acosh" 1)
  70 (def-math-rtn "atanh" 1)
  71
  72 ;;; Exponential and Logarithmic.
  73 #!-x86 (def-math-rtn "exp" 1)
  74 #!-x86 (def-math-rtn "log" 1)
  75 #!-x86 (def-math-rtn "log10" 1)
  76 (def-math-rtn "pow" 2)
  77 #!-x86 (def-math-rtn "sqrt" 1)
  78 (def-math-rtn "hypot" 2)
  79 #!-(or hpux x86) (def-math-rtn "log1p" 1)
  80
  81 #!+x86 ;; These are needed for use by byte-compiled files.
  82 (progn
  83   (defun %sin (x)
  84     (declare (double-float x)
  85              (values double-float))
  86     (%sin x))
  87   (defun %sin-quick (x)
  88     (declare (double-float x)
  89              (values double-float))
  90     (%sin-quick x))
  91   (defun %cos (x)
  92     (declare (double-float x)
  93              (values double-float))
  94     (%cos x))
  95   (defun %cos-quick (x)
  96     (declare (double-float x)
  97              (values double-float))
  98     (%cos-quick x))
  99   (defun %tan (x)
 100     (declare (double-float x)
 101              (values double-float))
 102     (%tan x))
 103   (defun %tan-quick (x)
 104     (declare (double-float x)
 105              (values double-float))
 106     (%tan-quick x))
 107   (defun %atan (x)
 108     (declare (double-float x)
 109              (values double-float))
 110     (%atan x))
 111   (defun %atan2 (x y)
 112     (declare (double-float x y)
 113              (values double-float))
 114     (%atan2 x y))
 115   (defun %exp (x)
 116     (declare (double-float x)
 117              (values double-float))
 118     (%exp x))
 119   (defun %log (x)
 120     (declare (double-float x)
 121              (values double-float))
 122     (%log x))
 123   (defun %log10 (x)
 124     (declare (double-float x)
 125              (values double-float))
 126     (%log10 x))
 127   #+nil ;; notyet
 128   (defun %pow (x y)
 129     (declare (type (double-float 0d0) x)
 130              (double-float y)
 131              (values (double-float 0d0)))
 132     (%pow x y))
 133   (defun %sqrt (x)
 134     (declare (double-float x)
 135              (values double-float))
 136     (%sqrt x))
 137   (defun %scalbn (f ex)
 138     (declare (double-float f)
 139              (type (signed-byte 32) ex)
 140              (values double-float))
 141     (%scalbn f ex))
 142   (defun %scalb (f ex)
 143     (declare (double-float f ex)
 144              (values double-float))
 145     (%scalb f ex))
 146   (defun %logb (x)
 147     (declare (double-float x)
 148              (values double-float))
 149     (%logb x))
 150   (defun %log1p (x)
 151     (declare (double-float x)
 152              (values double-float))
 153     (%log1p x))
 154   ) ; progn
 155 \f
 156 ;;;; power functions
 157
 158 (defun exp (number)
 159   #!+sb-doc
 160   "Return e raised to the power NUMBER."
 161   (number-dispatch ((number number))
 162     (handle-reals %exp number)
 163     ((complex)
 164      (* (exp (realpart number))
 165         (cis (imagpart number))))))
 166
 167 ;;; INTEXP -- Handle the rational base, integer power case.
 168
 169 ;;; FIXME: As long as the
 170 ;;; system dies on stack overflow or memory exhaustion, it seems reasonable
 171 ;;; to have this, but its default should be NIL, and when it's NIL,
 172 ;;; anything should be accepted.
 173 (defparameter *intexp-maximum-exponent* 10000)
 174
 175 ;;; This function precisely calculates base raised to an integral power. It
 176 ;;; separates the cases by the sign of power, for efficiency reasons, as powers
 177 ;;; can be calculated more efficiently if power is a positive integer. Values
 178 ;;; of power are calculated as positive integers, and inverted if negative.
 179 (defun intexp (base power)
 180   (when (> (abs power) *intexp-maximum-exponent*)
 181     ;; FIXME: should be ordinary error, not CERROR. (Once we set the
 182     ;; default for the variable to NIL, the un-continuable error will
 183     ;; be less obnoxious.)
 184     (cerror "Continue with calculation."
 185             "The absolute value of ~S exceeds ~S."
 186             power '*intexp-maximum-exponent* base power))
 187   (cond ((minusp power)
 188          (/ (intexp base (- power))))
 189         ((eql base 2)
 190          (ash 1 power))
 191         (t
 192          (do ((nextn (ash power -1) (ash power -1))
 193               (total (if (oddp power) base 1)
 194                      (if (oddp power) (* base total) total)))
 195              ((zerop nextn) total)
 196            (setq base (* base base))
 197            (setq power nextn)))))
 198
 199 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
 200 ;;; floating point stuff. If both args are real, we try %POW right off,
 201 ;;; assuming it will return 0 if the result may be complex. If so, we call
 202 ;;; COMPLEX-POW which directly computes the complex result. We also separate
 203 ;;; the complex-real and real-complex cases from the general complex case.
 204 (defun expt (base power)
 205   #!+sb-doc
 206   "Returns BASE raised to the POWER."
 207   (if (zerop power)
 208       (1+ (* base power))
 209     (labels (;; determine if the double float is an integer.
 210              ;;  0 - not an integer
 211              ;;  1 - an odd int
 212              ;;  2 - an even int
 213              (isint (ihi lo)
 214                (declare (type (unsigned-byte 31) ihi)
 215                         (type (unsigned-byte 32) lo)
 216                         (optimize (speed 3) (safety 0)))
 217                (let ((isint 0))
 218                  (declare (type fixnum isint))
 219                  (cond ((>= ihi #x43400000)     ; exponent >= 53
 220                         (setq isint 2))
 221                        ((>= ihi #x3ff00000)
 222                         (let ((k (- (ash ihi -20) #x3ff)))      ; exponent
 223                           (declare (type (mod 53) k))
 224                           (cond ((> k 20)
 225                                  (let* ((shift (- 52 k))
 226                                         (j (logand (ash lo (- shift))))
 227                                         (j2 (ash j shift)))
 228                                    (declare (type (mod 32) shift)
 229                                             (type (unsigned-byte 32) j j2))
 230                                    (when (= j2 lo)
 231                                      (setq isint (- 2 (logand j 1))))))
 232                                 ((= lo 0)
 233                                  (let* ((shift (- 20 k))
 234                                         (j (ash ihi (- shift)))
 235                                         (j2 (ash j shift)))
 236                                    (declare (type (mod 32) shift)
 237                                             (type (unsigned-byte 31) j j2))
 238                                    (when (= j2 ihi)
 239                                      (setq isint (- 2 (logand j 1))))))))))
 240                  isint))
 241              (real-expt (x y rtype)
 242                (let ((x (coerce x 'double-float))
 243                      (y (coerce y 'double-float)))
 244                  (declare (double-float x y))
 245                  (let* ((x-hi (sb!kernel:double-float-high-bits x))
 246                         (x-lo (sb!kernel:double-float-low-bits x))
 247                         (x-ihi (logand x-hi #x7fffffff))
 248                         (y-hi (sb!kernel:double-float-high-bits y))
 249                         (y-lo (sb!kernel:double-float-low-bits y))
 250                         (y-ihi (logand y-hi #x7fffffff)))
 251                    (declare (type (signed-byte 32) x-hi y-hi)
 252                             (type (unsigned-byte 31) x-ihi y-ihi)
 253                             (type (unsigned-byte 32) x-lo y-lo))
 254                    ;; y==zero: x**0 = 1
 255                    (when (zerop (logior y-ihi y-lo))
 256                      (return-from real-expt (coerce 1d0 rtype)))
 257                    ;; +-NaN return x+y
 258                    (when (or (> x-ihi #x7ff00000)
 259                              (and (= x-ihi #x7ff00000) (/= x-lo 0))
 260                              (> y-ihi #x7ff00000)
 261                              (and (= y-ihi #x7ff00000) (/= y-lo 0)))
 262                      (return-from real-expt (coerce (+ x y) rtype)))
 263                    (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
 264                      (declare (type fixnum yisint))
 265                      ;; special value of y
 266                      (when (and (zerop y-lo) (= y-ihi #x7ff00000))
 267                        ;; y is +-inf
 268                        (return-from real-expt
 269                          (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
 270                                 ;; +-1**inf is NaN
 271                                 (coerce (- y y) rtype))
 272                                ((>= x-ihi #x3ff00000)
 273                                 ;; (|x|>1)**+-inf = inf,0
 274                                 (if (>= y-hi 0)
 275                                     (coerce y rtype)
 276                                     (coerce 0 rtype)))
 277                                (t
 278                                 ;; (|x|<1)**-,+inf = inf,0
 279                                 (if (< y-hi 0)
 280                                     (coerce (- y) rtype)
 281                                     (coerce 0 rtype))))))
 282
 283                      (let ((abs-x (abs x)))
 284                        (declare (double-float abs-x))
 285                        ;; special value of x
 286                        (when (and (zerop x-lo)
 287                                   (or (= x-ihi #x7ff00000) (zerop x-ihi)
 288                                       (= x-ihi #x3ff00000)))
 289                          ;; x is +-0,+-inf,+-1
 290                          (let ((z (if (< y-hi 0)
 291                                       (/ 1 abs-x)       ; z = (1/|x|)
 292                                       abs-x)))
 293                            (declare (double-float z))
 294                            (when (< x-hi 0)
 295                              (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
 296                                     ;; (-1)**non-int
 297                                     (let ((y*pi (* y pi)))
 298                                       (declare (double-float y*pi))
 299                                       (return-from real-expt
 300                                         (complex
 301                                          (coerce (%cos y*pi) rtype)
 302                                          (coerce (%sin y*pi) rtype)))))
 303                                    ((= yisint 1)
 304                                     ;; (x<0)**odd = -(|x|**odd)
 305                                     (setq z (- z)))))
 306                            (return-from real-expt (coerce z rtype))))
 307
 308                        (if (>= x-hi 0)
 309                            ;; x>0
 310                            (coerce (sb!kernel::%pow x y) rtype)
 311                            ;; x<0
 312                            (let ((pow (sb!kernel::%pow abs-x y)))
 313                              (declare (double-float pow))
 314                              (case yisint
 315                                (1 ; Odd
 316                                 (coerce (* -1d0 pow) rtype))
 317                                (2 ; Even
 318                                 (coerce pow rtype))
 319                                (t ; Non-integer
 320                                 (let ((y*pi (* y pi)))
 321                                   (declare (double-float y*pi))
 322                                   (complex
 323                                    (coerce (* pow (%cos y*pi)) rtype)
 324                                    (coerce (* pow (%sin y*pi)) rtype)))))))))))))
 325       (declare (inline real-expt))
 326       (number-dispatch ((base number) (power number))
 327         (((foreach fixnum (or bignum ratio) (complex rational)) integer)
 328          (intexp base power))
 329         (((foreach single-float double-float) rational)
 330          (real-expt base power '(dispatch-type base)))
 331         (((foreach fixnum (or bignum ratio) single-float)
 332           (foreach ratio single-float))
 333          (real-expt base power 'single-float))
 334         (((foreach fixnum (or bignum ratio) single-float double-float)
 335           double-float)
 336          (real-expt base power 'double-float))
 337         ((double-float single-float)
 338          (real-expt base power 'double-float))
 339         (((foreach (complex rational) (complex float)) rational)
 340          (* (expt (abs base) power)
 341             (cis (* power (phase base)))))
 342         (((foreach fixnum (or bignum ratio) single-float double-float)
 343           complex)
 344          (if (and (zerop base) (plusp (realpart power)))
 345              (* base power)
 346              (exp (* power (log base)))))
 347         (((foreach (complex float) (complex rational))
 348           (foreach complex double-float single-float))
 349          (if (and (zerop base) (plusp (realpart power)))
 350              (* base power)
 351              (exp (* power (log base)))))))))
 352
 353 (defun log (number &optional (base nil base-p))
 354   #!+sb-doc
 355   "Return the logarithm of NUMBER in the base BASE, which defaults to e."
 356   (if base-p
 357       (if (zerop base)
 358           base                          ; ANSI spec
 359           (/ (log number) (log base)))
 360       (number-dispatch ((number number))
 361         (((foreach fixnum bignum ratio))
 362          (if (minusp number)
 363              (complex (log (- number)) (coerce pi 'single-float))
 364              (coerce (%log (coerce number 'double-float)) 'single-float)))
 365         (((foreach single-float double-float))
 366          ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
 367          ;; Since this doesn't seem to be an implementation issue
 368          ;; I (pw) take the Kahan result.
 369          (if (< (float-sign number)
 370                 (coerce 0 '(dispatch-type number)))
 371              (complex (log (- number)) (coerce pi '(dispatch-type number)))
 372              (coerce (%log (coerce number 'double-float))
 373                      '(dispatch-type number))))
 374         ((complex)
 375          (complex-log number)))))
 376
 377 (defun sqrt (number)
 378   #!+sb-doc
 379   "Return the square root of NUMBER."
 380   (number-dispatch ((number number))
 381     (((foreach fixnum bignum ratio))
 382      (if (minusp number)
 383          (complex-sqrt number)
 384          (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
 385     (((foreach single-float double-float))
 386      (if (minusp number)
 387          (complex-sqrt number)
 388          (coerce (%sqrt (coerce number 'double-float))
 389                  '(dispatch-type number))))
 390      ((complex)
 391       (complex-sqrt number))))
 392 \f
 393 ;;;; trigonometic and related functions
 394
 395 (defun abs (number)
 396   #!+sb-doc
 397   "Returns the absolute value of the number."
 398   (number-dispatch ((number number))
 399     (((foreach single-float double-float fixnum rational))
 400      (abs number))
 401     ((complex)
 402      (let ((rx (realpart number))
 403            (ix (imagpart number)))
 404        (etypecase rx
 405          (rational
 406           (sqrt (+ (* rx rx) (* ix ix))))
 407          (single-float
 408           (coerce (%hypot (coerce rx 'double-float)
 409                           (coerce ix 'double-float))
 410                   'single-float))
 411          (double-float
 412           (%hypot rx ix)))))))
 413
 414 (defun phase (number)
 415   #!+sb-doc
 416   "Returns the angle part of the polar representation of a complex number.
 417   For complex numbers, this is (atan (imagpart number) (realpart number)).
 418   For non-complex positive numbers, this is 0. For non-complex negative
 419   numbers this is PI."
 420   (etypecase number
 421     (rational
 422      (if (minusp number)
 423          (coerce pi 'single-float)
 424          0.0f0))
 425     (single-float
 426      (if (minusp (float-sign number))
 427          (coerce pi 'single-float)
 428          0.0f0))
 429     (double-float
 430      (if (minusp (float-sign number))
 431          (coerce pi 'double-float)
 432          0.0d0))
 433     (complex
 434      (atan (imagpart number) (realpart number)))))
 435
 436 (defun sin (number)
 437   #!+sb-doc
 438   "Return the sine of NUMBER."
 439   (number-dispatch ((number number))
 440     (handle-reals %sin number)
 441     ((complex)
 442      (let ((x (realpart number))
 443            (y (imagpart number)))
 444        (complex (* (sin x) (cosh y))
 445                 (* (cos x) (sinh y)))))))
 446
 447 (defun cos (number)
 448   #!+sb-doc
 449   "Return the cosine of NUMBER."
 450   (number-dispatch ((number number))
 451     (handle-reals %cos number)
 452     ((complex)
 453      (let ((x (realpart number))
 454            (y (imagpart number)))
 455        (complex (* (cos x) (cosh y))
 456                 (- (* (sin x) (sinh y))))))))
 457
 458 (defun tan (number)
 459   #!+sb-doc
 460   "Return the tangent of NUMBER."
 461   (number-dispatch ((number number))
 462     (handle-reals %tan number)
 463     ((complex)
 464      (complex-tan number))))
 465
 466 (defun cis (theta)
 467   #!+sb-doc
 468   "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
 469   (if (complexp theta)
 470       (error "Argument to CIS is complex: ~S" theta)
 471       (complex (cos theta) (sin theta))))
 472
 473 (defun asin (number)
 474   #!+sb-doc
 475   "Return the arc sine of NUMBER."
 476   (number-dispatch ((number number))
 477     ((rational)
 478      (if (or (> number 1) (< number -1))
 479          (complex-asin number)
 480          (coerce (%asin (coerce number 'double-float)) 'single-float)))
 481     (((foreach single-float double-float))
 482      (if (or (> number (coerce 1 '(dispatch-type number)))
 483              (< number (coerce -1 '(dispatch-type number))))
 484          (complex-asin number)
 485          (coerce (%asin (coerce number 'double-float))
 486                  '(dispatch-type number))))
 487     ((complex)
 488      (complex-asin number))))
 489
 490 (defun acos (number)
 491   #!+sb-doc
 492   "Return the arc cosine of NUMBER."
 493   (number-dispatch ((number number))
 494     ((rational)
 495      (if (or (> number 1) (< number -1))
 496          (complex-acos number)
 497          (coerce (%acos (coerce number 'double-float)) 'single-float)))
 498     (((foreach single-float double-float))
 499      (if (or (> number (coerce 1 '(dispatch-type number)))
 500              (< number (coerce -1 '(dispatch-type number))))
 501          (complex-acos number)
 502          (coerce (%acos (coerce number 'double-float))
 503                  '(dispatch-type number))))
 504     ((complex)
 505      (complex-acos number))))
 506
 507 (defun atan (y &optional (x nil xp))
 508   #!+sb-doc
 509   "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
 510   (if xp
 511       (flet ((atan2 (y x)
 512                (declare (type double-float y x)
 513                         (values double-float))
 514                (if (zerop x)
 515                    (if (zerop y)
 516                        (if (plusp (float-sign x))
 517                            y
 518                            (float-sign y pi))
 519                        (float-sign y (/ pi 2)))
 520                    (%atan2 y x))))
 521         (number-dispatch ((y number) (x number))
 522           ((double-float
 523             (foreach double-float single-float fixnum bignum ratio))
 524            (atan2 y (coerce x 'double-float)))
 525           (((foreach single-float fixnum bignum ratio)
 526             double-float)
 527            (atan2 (coerce y 'double-float) x))
 528           (((foreach single-float fixnum bignum ratio)
 529             (foreach single-float fixnum bignum ratio))
 530            (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
 531                    'single-float))))
 532       (number-dispatch ((y number))
 533         (handle-reals %atan y)
 534         ((complex)
 535          (complex-atan y)))))
 536
 537 ;; It seems that everyone has a C version of sinh, cosh, and
 538 ;; tanh. Let's use these for reals because the original
 539 ;; implementations based on the definitions lose big in round-off
 540 ;; error. These bad definitions also mean that sin and cos for
 541 ;; complex numbers can also lose big.
 542
 543 #+nil
 544 (defun sinh (number)
 545   #!+sb-doc
 546   "Return the hyperbolic sine of NUMBER."
 547   (/ (- (exp number) (exp (- number))) 2))
 548
 549 (defun sinh (number)
 550   #!+sb-doc
 551   "Return the hyperbolic sine of NUMBER."
 552   (number-dispatch ((number number))
 553     (handle-reals %sinh number)
 554     ((complex)
 555      (let ((x (realpart number))
 556            (y (imagpart number)))
 557        (complex (* (sinh x) (cos y))
 558                 (* (cosh x) (sin y)))))))
 559
 560 #+nil
 561 (defun cosh (number)
 562   #!+sb-doc
 563   "Return the hyperbolic cosine of NUMBER."
 564   (/ (+ (exp number) (exp (- number))) 2))
 565
 566 (defun cosh (number)
 567   #!+sb-doc
 568   "Return the hyperbolic cosine of NUMBER."
 569   (number-dispatch ((number number))
 570     (handle-reals %cosh number)
 571     ((complex)
 572      (let ((x (realpart number))
 573            (y (imagpart number)))
 574        (complex (* (cosh x) (cos y))
 575                 (* (sinh x) (sin y)))))))
 576
 577 (defun tanh (number)
 578   #!+sb-doc
 579   "Return the hyperbolic tangent of NUMBER."
 580   (number-dispatch ((number number))
 581     (handle-reals %tanh number)
 582     ((complex)
 583      (complex-tanh number))))
 584
 585 (defun asinh (number)
 586   #!+sb-doc
 587   "Return the hyperbolic arc sine of NUMBER."
 588   (number-dispatch ((number number))
 589     (handle-reals %asinh number)
 590     ((complex)
 591      (complex-asinh number))))
 592
 593 (defun acosh (number)
 594   #!+sb-doc
 595   "Return the hyperbolic arc cosine of NUMBER."
 596   (number-dispatch ((number number))
 597     ((rational)
 598      ;; acosh is complex if number < 1
 599      (if (< number 1)
 600          (complex-acosh number)
 601          (coerce (%acosh (coerce number 'double-float)) 'single-float)))
 602     (((foreach single-float double-float))
 603      (if (< number (coerce 1 '(dispatch-type number)))
 604          (complex-acosh number)
 605          (coerce (%acosh (coerce number 'double-float))
 606                  '(dispatch-type number))))
 607     ((complex)
 608      (complex-acosh number))))
 609
 610 (defun atanh (number)
 611   #!+sb-doc
 612   "Return the hyperbolic arc tangent of NUMBER."
 613   (number-dispatch ((number number))
 614     ((rational)
 615      ;; atanh is complex if |number| > 1
 616      (if (or (> number 1) (< number -1))
 617          (complex-atanh number)
 618          (coerce (%atanh (coerce number 'double-float)) 'single-float)))
 619     (((foreach single-float double-float))
 620      (if (or (> number (coerce 1 '(dispatch-type number)))
 621              (< number (coerce -1 '(dispatch-type number))))
 622          (complex-atanh number)
 623          (coerce (%atanh (coerce number 'double-float))
 624                  '(dispatch-type number))))
 625     ((complex)
 626      (complex-atanh number))))
 627
 628 ;;; HP-UX does not supply a C version of log1p, so
 629 ;;; use the definition.
 630
 631 #!+hpux
 632 #!-sb-fluid (declaim (inline %log1p))
 633 #!+hpux
 634 (defun %log1p (number)
 635   (declare (double-float number)
 636            (optimize (speed 3) (safety 0)))
 637   (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
 638