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