src/compiler/float-tran.lisp

   1 ;;;; This file contains floating-point-specific transforms, and may be
   2 ;;;; somewhat implementation-dependent in its assumptions of what the
   3 ;;;; formats are.
   4
   5 ;;;; This software is part of the SBCL system. See the README file for
   6 ;;;; more information.
   7 ;;;;
   8 ;;;; This software is derived from the CMU CL system, which was
   9 ;;;; written at Carnegie Mellon University and released into the
  10 ;;;; public domain. The software is in the public domain and is
  11 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
  12 ;;;; files for more information.
  13
  14 (in-package "SB!C")
  15 \f
  16 ;;;; coercions
  17
  18 (defknown %single-float (real) single-float
  19     (movable foldable))
  20 (defknown %double-float (real) double-float
  21     (movable foldable))
  22
  23 (deftransform float ((n f) (* single-float) *)
  24   '(%single-float n))
  25
  26 (deftransform float ((n f) (* double-float) *)
  27   '(%double-float n))
  28
  29 (deftransform float ((n) *)
  30   '(if (floatp n)
  31        n
  32        (%single-float n)))
  33
  34 (deftransform %single-float ((n) (single-float) *)
  35   'n)
  36
  37 (deftransform %double-float ((n) (double-float) *)
  38   'n)
  39
  40 ;;; RANDOM
  41 (macrolet ((frob (fun type)
  42              `(deftransform random ((num &optional state)
  43                                     (,type &optional *) *)
  44                 "Use inline float operations."
  45                 '(,fun num (or state *random-state*)))))
  46   (frob %random-single-float single-float)
  47   (frob %random-double-float double-float))
  48
  49 ;;; Return an expression to generate an integer of N-BITS many random
  50 ;;; bits, using the minimal number of random chunks possible.
  51 (defun generate-random-expr-for-power-of-2 (n-bits state)
  52   (declare (type (integer 1 #.sb!vm:n-word-bits) n-bits))
  53   (multiple-value-bind (n-chunk-bits chunk-expr)
  54       (cond ((<= n-bits n-random-chunk-bits)
  55              (values n-random-chunk-bits `(random-chunk ,state)))
  56             ((<= n-bits (* 2 n-random-chunk-bits))
  57              (values (* 2 n-random-chunk-bits) `(big-random-chunk ,state)))
  58             (t
  59              (error "Unexpectedly small N-RANDOM-CHUNK-BITS")))
  60     (if (< n-bits n-chunk-bits)
  61         `(logand ,(1- (ash 1 n-bits)) ,chunk-expr)
  62         chunk-expr)))
  63
  64 ;;; This transform for compile-time constant word-sized integers
  65 ;;; generates an accept-reject loop to achieve equidistribution of the
  66 ;;; returned values. Several optimizations are done: If NUM is a power
  67 ;;; of two no loop is needed. If the random chunk size is half the word
  68 ;;; size only one chunk is used where sufficient. For values of NUM
  69 ;;; where it is possible and results in faster code, the rejection
  70 ;;; probability is reduced by accepting all values below the largest
  71 ;;; multiple of the limit that fits into one or two chunks and and doing
  72 ;;; a division to get the random value into the desired range.
  73 (deftransform random ((num &optional state)
  74                       ((constant-arg (integer 1 #.(expt 2 sb!vm:n-word-bits)))
  75                        &optional *)
  76                       *
  77                       :policy (and (> speed compilation-speed)
  78                                    (> speed space)))
  79   "optimize to inlined RANDOM-CHUNK operations"
  80   (let ((num (lvar-value num)))
  81     (if (= num 1)
  82         0
  83         (flet ((chunk-n-bits-and-expr (n-bits)
  84                  (cond ((<= n-bits n-random-chunk-bits)
  85                         (values n-random-chunk-bits
  86                                 '(random-chunk (or state *random-state*))))
  87                        ((<= n-bits (* 2 n-random-chunk-bits))
  88                         (values (* 2 n-random-chunk-bits)
  89                                 '(big-random-chunk (or state *random-state*))))
  90                        (t
  91                         (error "Unexpectedly small N-RANDOM-CHUNK-BITS")))))
  92           (if (zerop (logand num (1- num)))
  93               ;; NUM is a power of 2.
  94               (let ((n-bits (integer-length (1- num))))
  95                 (multiple-value-bind (n-chunk-bits chunk-expr)
  96                     (chunk-n-bits-and-expr n-bits)
  97                   (if (< n-bits n-chunk-bits)
  98                       `(logand ,(1- (ash 1 n-bits)) ,chunk-expr)
  99                       chunk-expr)))
 100               ;; Generate an accept-reject loop.
 101               (let ((n-bits (integer-length num)))
 102                 (multiple-value-bind (n-chunk-bits chunk-expr)
 103                     (chunk-n-bits-and-expr n-bits)
 104                   (if (or (> (* num 3) (expt 2 n-chunk-bits))
 105                           (logbitp (- n-bits 2) num))
 106                       ;; Division can't help as the quotient is below 3,
 107                       ;; or is too costly as the rejection probability
 108                       ;; without it is already small (namely at most 1/4
 109                       ;; with the given test, which is experimentally a
 110                       ;; reasonable threshold and cheap to test for).
 111                       `(loop
 112                          (let ((bits ,(generate-random-expr-for-power-of-2
 113                                        n-bits '(or state *random-state*))))
 114                            (when (< bits num)
 115                              (return bits))))
 116                       (let ((d (truncate (expt 2 n-chunk-bits) num)))
 117                         `(loop
 118                            (let ((bits ,chunk-expr))
 119                              (when (< bits ,(* num d))
 120                                (return (values (truncate bits ,d)))))))))))))))
 121
 122 \f
 123 ;;;; float accessors
 124
 125 (defknown make-single-float ((signed-byte 32)) single-float
 126   (movable flushable))
 127
 128 (defknown make-double-float ((signed-byte 32) (unsigned-byte 32)) double-float
 129   (movable flushable))
 130
 131 #-sb-xc-host
 132 (deftransform make-single-float ((bits)
 133                                  ((signed-byte 32)))
 134   "Conditional constant folding"
 135   (unless (constant-lvar-p bits)
 136     (give-up-ir1-transform))
 137   (let* ((bits  (lvar-value bits))
 138          (float (make-single-float bits)))
 139     (when (float-nan-p float)
 140       (give-up-ir1-transform))
 141     float))
 142
 143 #-sb-xc-host
 144 (deftransform make-double-float ((hi lo)
 145                                  ((signed-byte 32) (unsigned-byte 32)))
 146   "Conditional constant folding"
 147   (unless (and (constant-lvar-p hi)
 148                (constant-lvar-p lo))
 149     (give-up-ir1-transform))
 150   (let* ((hi    (lvar-value hi))
 151          (lo    (lvar-value lo))
 152          (float (make-double-float hi lo)))
 153     (when (float-nan-p float)
 154       (give-up-ir1-transform))
 155     float))
 156
 157 (defknown single-float-bits (single-float) (signed-byte 32)
 158   (movable foldable flushable))
 159
 160 (defknown double-float-high-bits (double-float) (signed-byte 32)
 161   (movable foldable flushable))
 162
 163 (defknown double-float-low-bits (double-float) (unsigned-byte 32)
 164   (movable foldable flushable))
 165
 166 (deftransform float-sign ((float &optional float2)
 167                           (single-float &optional single-float) *)
 168   (if float2
 169       (let ((temp (gensym)))
 170         `(let ((,temp (abs float2)))
 171           (if (minusp (single-float-bits float)) (- ,temp) ,temp)))
 172       '(if (minusp (single-float-bits float)) -1f0 1f0)))
 173
 174 (deftransform float-sign ((float &optional float2)
 175                           (double-float &optional double-float) *)
 176   (if float2
 177       (let ((temp (gensym)))
 178         `(let ((,temp (abs float2)))
 179           (if (minusp (double-float-high-bits float)) (- ,temp) ,temp)))
 180       '(if (minusp (double-float-high-bits float)) -1d0 1d0)))
 181 \f
 182 ;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, and SCALE-FLOAT
 183
 184 (defknown decode-single-float (single-float)
 185   (values single-float single-float-exponent (single-float -1f0 1f0))
 186   (movable foldable flushable))
 187
 188 (defknown decode-double-float (double-float)
 189   (values double-float double-float-exponent (double-float -1d0 1d0))
 190   (movable foldable flushable))
 191
 192 (defknown integer-decode-single-float (single-float)
 193   (values single-float-significand single-float-int-exponent (integer -1 1))
 194   (movable foldable flushable))
 195
 196 (defknown integer-decode-double-float (double-float)
 197   (values double-float-significand double-float-int-exponent (integer -1 1))
 198   (movable foldable flushable))
 199
 200 (defknown scale-single-float (single-float integer) single-float
 201   (movable foldable flushable))
 202
 203 (defknown scale-double-float (double-float integer) double-float
 204   (movable foldable flushable))
 205
 206 (deftransform decode-float ((x) (single-float) *)
 207   '(decode-single-float x))
 208
 209 (deftransform decode-float ((x) (double-float) *)
 210   '(decode-double-float x))
 211
 212 (deftransform integer-decode-float ((x) (single-float) *)
 213   '(integer-decode-single-float x))
 214
 215 (deftransform integer-decode-float ((x) (double-float) *)
 216   '(integer-decode-double-float x))
 217
 218 (deftransform scale-float ((f ex) (single-float *) *)
 219   (if (and #!+x86 t #!-x86 nil
 220            (csubtypep (lvar-type ex)
 221                       (specifier-type '(signed-byte 32))))
 222       '(coerce (%scalbn (coerce f 'double-float) ex) 'single-float)
 223       '(scale-single-float f ex)))
 224
 225 (deftransform scale-float ((f ex) (double-float *) *)
 226   (if (and #!+x86 t #!-x86 nil
 227            (csubtypep (lvar-type ex)
 228                       (specifier-type '(signed-byte 32))))
 229       '(%scalbn f ex)
 230       '(scale-double-float f ex)))
 231
 232 ;;; What is the CROSS-FLOAT-INFINITY-KLUDGE?
 233 ;;;
 234 ;;; SBCL's own implementation of floating point supports floating
 235 ;;; point infinities. Some of the old CMU CL :PROPAGATE-FLOAT-TYPE and
 236 ;;; :PROPAGATE-FUN-TYPE code, like the DEFOPTIMIZERs below, uses this
 237 ;;; floating point support. Thus, we have to avoid running it on the
 238 ;;; cross-compilation host, since we're not guaranteed that the
 239 ;;; cross-compilation host will support floating point infinities.
 240 ;;;
 241 ;;; If we wanted to live dangerously, we could conditionalize the code
 242 ;;; with #+(OR SBCL SB-XC) instead. That way, if the cross-compilation
 243 ;;; host happened to be SBCL, we'd be able to run the infinity-using
 244 ;;; code. Pro:
 245 ;;;   * SBCL itself gets built with more complete optimization.
 246 ;;; Con:
 247 ;;;   * You get a different SBCL depending on what your cross-compilation
 248 ;;;     host is.
 249 ;;; So far the pros and cons seem seem to be mostly academic, since
 250 ;;; AFAIK (WHN 2001-08-28) the propagate-foo-type optimizations aren't
 251 ;;; actually important in compiling SBCL itself. If this changes, then
 252 ;;; we have to decide:
 253 ;;;   * Go for simplicity, leaving things as they are.
 254 ;;;   * Go for performance at the expense of conceptual clarity,
 255 ;;;     using #+(OR SBCL SB-XC) and otherwise leaving the build
 256 ;;;     process as is.
 257 ;;;   * Go for performance at the expense of build time, using
 258 ;;;     #+(OR SBCL SB-XC) and also making SBCL do not just
 259 ;;;     make-host-1.sh and make-host-2.sh, but a third step
 260 ;;;     make-host-3.sh where it builds itself under itself. (Such a
 261 ;;;     3-step build process could also help with other things, e.g.
 262 ;;;     using specialized arrays to represent debug information.)
 263 ;;;   * Rewrite the code so that it doesn't depend on unportable
 264 ;;;     floating point infinities.
 265
 266 ;;; optimizers for SCALE-FLOAT. If the float has bounds, new bounds
 267 ;;; are computed for the result, if possible.
 268 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 269 (progn
 270
 271 (defun scale-float-derive-type-aux (f ex same-arg)
 272   (declare (ignore same-arg))
 273   (flet ((scale-bound (x n)
 274            ;; We need to be a bit careful here and catch any overflows
 275            ;; that might occur. We can ignore underflows which become
 276            ;; zeros.
 277            (set-bound
 278             (handler-case
 279              (scale-float (type-bound-number x) n)
 280              (floating-point-overflow ()
 281                 nil))
 282             (consp x))))
 283     (when (and (numeric-type-p f) (numeric-type-p ex))
 284       (let ((f-lo (numeric-type-low f))
 285             (f-hi (numeric-type-high f))
 286             (ex-lo (numeric-type-low ex))
 287             (ex-hi (numeric-type-high ex))
 288             (new-lo nil)
 289             (new-hi nil))
 290         (when f-hi
 291           (if (< (float-sign (type-bound-number f-hi)) 0.0)
 292               (when ex-lo
 293                 (setf new-hi (scale-bound f-hi ex-lo)))
 294               (when ex-hi
 295                 (setf new-hi (scale-bound f-hi ex-hi)))))
 296         (when f-lo
 297           (if (< (float-sign (type-bound-number f-lo)) 0.0)
 298               (when ex-hi
 299                 (setf new-lo (scale-bound f-lo ex-hi)))
 300               (when ex-lo
 301                 (setf new-lo (scale-bound f-lo ex-lo)))))
 302         (make-numeric-type :class (numeric-type-class f)
 303                            :format (numeric-type-format f)
 304                            :complexp :real
 305                            :low new-lo
 306                            :high new-hi)))))
 307 (defoptimizer (scale-single-float derive-type) ((f ex))
 308   (two-arg-derive-type f ex #'scale-float-derive-type-aux
 309                        #'scale-single-float t))
 310 (defoptimizer (scale-double-float derive-type) ((f ex))
 311   (two-arg-derive-type f ex #'scale-float-derive-type-aux
 312                        #'scale-double-float t))
 313
 314 ;;; DEFOPTIMIZERs for %SINGLE-FLOAT and %DOUBLE-FLOAT. This makes the
 315 ;;; FLOAT function return the correct ranges if the input has some
 316 ;;; defined range. Quite useful if we want to convert some type of
 317 ;;; bounded integer into a float.
 318 (macrolet
 319     ((frob (fun type most-negative most-positive)
 320        (let ((aux-name (symbolicate fun "-DERIVE-TYPE-AUX")))
 321          `(progn
 322             (defun ,aux-name (num)
 323               ;; When converting a number to a float, the limits are
 324               ;; the same.
 325               (let* ((lo (bound-func (lambda (x)
 326                                        (if (< x ,most-negative)
 327                                            ,most-negative
 328                                            (coerce x ',type)))
 329                                      (numeric-type-low num)))
 330                      (hi (bound-func (lambda (x)
 331                                        (if (< ,most-positive x )
 332                                            ,most-positive
 333                                            (coerce x ',type)))
 334                                      (numeric-type-high num))))
 335                 (specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
 336
 337             (defoptimizer (,fun derive-type) ((num))
 338               (handler-case
 339                   (one-arg-derive-type num #',aux-name #',fun)
 340                 (type-error ()
 341                   nil)))))))
 342   (frob %single-float single-float
 343         most-negative-single-float most-positive-single-float)
 344   (frob %double-float double-float
 345         most-negative-double-float most-positive-double-float))
 346 ) ; PROGN
 347 \f
 348 ;;;; float contagion
 349
 350 (defun safe-ctype-for-single-coercion-p (x)
 351   ;; See comment in SAFE-SINGLE-COERCION-P -- this deals with the same
 352   ;; problem, but in the context of evaluated and compiled (+ <int> <single>)
 353   ;; giving different result if we fail to check for this.
 354   (or (not (csubtypep x (specifier-type 'integer)))
 355       #!+x86
 356       (csubtypep x (specifier-type `(integer ,most-negative-exactly-single-float-fixnum
 357                                              ,most-positive-exactly-single-float-fixnum)))
 358       #!-x86
 359       (csubtypep x (specifier-type 'fixnum))))
 360
 361 ;;; Do some stuff to recognize when the loser is doing mixed float and
 362 ;;; rational arithmetic, or different float types, and fix it up. If
 363 ;;; we don't, he won't even get so much as an efficiency note.
 364 (deftransform float-contagion-arg1 ((x y) * * :defun-only t :node node)
 365   (if (or (not (types-equal-or-intersect (lvar-type y) (specifier-type 'single-float)))
 366           (safe-ctype-for-single-coercion-p (lvar-type x)))
 367       `(,(lvar-fun-name (basic-combination-fun node))
 368          (float x y) y)
 369       (give-up-ir1-transform)))
 370 (deftransform float-contagion-arg2 ((x y) * * :defun-only t :node node)
 371   (if (or (not (types-equal-or-intersect (lvar-type x) (specifier-type 'single-float)))
 372           (safe-ctype-for-single-coercion-p (lvar-type y)))
 373       `(,(lvar-fun-name (basic-combination-fun node))
 374          x (float y x))
 375       (give-up-ir1-transform)))
 376
 377 (dolist (x '(+ * / -))
 378   (%deftransform x '(function (rational float) *) #'float-contagion-arg1)
 379   (%deftransform x '(function (float rational) *) #'float-contagion-arg2))
 380
 381 (dolist (x '(= < > + * / -))
 382   (%deftransform x '(function (single-float double-float) *)
 383                  #'float-contagion-arg1)
 384   (%deftransform x '(function (double-float single-float) *)
 385                  #'float-contagion-arg2))
 386
 387 (macrolet ((def (type &rest args)
 388              `(deftransform * ((x y) (,type (constant-arg (member ,@args))) *
 389                                ;; Beware the SNaN!
 390                                :policy (zerop float-accuracy))
 391                 "optimize multiplication by one"
 392                 (let ((y (lvar-value y)))
 393                   (if (minusp y)
 394                       '(%negate x)
 395                       'x)))))
 396   (def single-float 1.0 -1.0)
 397   (def double-float 1.0d0 -1.0d0))
 398
 399 ;;; Return the reciprocal of X if it can be represented exactly, NIL otherwise.
 400 (defun maybe-exact-reciprocal (x)
 401   (unless (zerop x)
 402     (handler-case
 403         (multiple-value-bind (significand exponent sign)
 404             (integer-decode-float x)
 405           ;; only powers of 2 can be inverted exactly
 406           (unless (zerop (logand significand (1- significand)))
 407             (return-from maybe-exact-reciprocal nil))
 408           (let ((expected   (/ sign significand (expt 2 exponent)))
 409                 (reciprocal (/ x)))
 410             (multiple-value-bind (significand exponent sign)
 411                 (integer-decode-float reciprocal)
 412               ;; Denorms can't be inverted safely.
 413               (and (eql expected (* sign significand (expt 2 exponent)))
 414                    reciprocal))))
 415       (error () (return-from maybe-exact-reciprocal nil)))))
 416
 417 ;;; Replace constant division by multiplication with exact reciprocal,
 418 ;;; if one exists.
 419 (macrolet ((def (type)
 420              `(deftransform / ((x y) (,type (constant-arg ,type)) *
 421                                :node node)
 422                 "convert to multiplication by reciprocal"
 423                 (let ((n (lvar-value y)))
 424                   (if (policy node (zerop float-accuracy))
 425                       `(* x ,(/ n))
 426                       (let ((r (maybe-exact-reciprocal n)))
 427                         (if r
 428                             `(* x ,r)
 429                             (give-up-ir1-transform
 430                              "~S does not have an exact reciprocal"
 431                              n))))))))
 432   (def single-float)
 433   (def double-float))
 434
 435 ;;; Optimize addition and subtraction of zero
 436 (macrolet ((def (op type &rest args)
 437              `(deftransform ,op ((x y) (,type (constant-arg (member ,@args))) *
 438                                  ;; Beware the SNaN!
 439                                  :policy (zerop float-accuracy))
 440                 'x)))
 441   ;; No signed zeros, thanks.
 442   (def + single-float 0 0.0)
 443   (def - single-float 0 0.0)
 444   (def + double-float 0 0.0 0.0d0)
 445   (def - double-float 0 0.0 0.0d0))
 446
 447 ;;; On most platforms (+ x x) is faster than (* x 2)
 448 (macrolet ((def (type &rest args)
 449              `(deftransform * ((x y) (,type (constant-arg (member ,@args))))
 450                 '(+ x x))))
 451   (def single-float 2 2.0)
 452   (def double-float 2 2.0 2.0d0))
 453
 454 ;;; Prevent ZEROP, PLUSP, and MINUSP from losing horribly. We can't in
 455 ;;; general float rational args to comparison, since Common Lisp
 456 ;;; semantics says we are supposed to compare as rationals, but we can
 457 ;;; do it for any rational that has a precise representation as a
 458 ;;; float (such as 0).
 459 (macrolet ((frob (op)
 460              `(deftransform ,op ((x y) (float rational) *)
 461                 "open-code FLOAT to RATIONAL comparison"
 462                 (unless (constant-lvar-p y)
 463                   (give-up-ir1-transform
 464                    "The RATIONAL value isn't known at compile time."))
 465                 (let ((val (lvar-value y)))
 466                   (unless (eql (rational (float val)) val)
 467                     (give-up-ir1-transform
 468                      "~S doesn't have a precise float representation."
 469                      val)))
 470                 `(,',op x (float y x)))))
 471   (frob <)
 472   (frob >)
 473   (frob =))
 474 \f
 475 ;;;; irrational derive-type methods
 476
 477 ;;; Derive the result to be float for argument types in the
 478 ;;; appropriate domain.
 479 #+sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 480 (dolist (stuff '((asin (real -1.0 1.0))
 481                  (acos (real -1.0 1.0))
 482                  (acosh (real 1.0))
 483                  (atanh (real -1.0 1.0))
 484                  (sqrt (real 0.0))))
 485   (destructuring-bind (name type) stuff
 486     (let ((type (specifier-type type)))
 487       (setf (fun-info-derive-type (fun-info-or-lose name))
 488             (lambda (call)
 489               (declare (type combination call))
 490               (when (csubtypep (lvar-type
 491                                 (first (combination-args call)))
 492                                type)
 493                 (specifier-type 'float)))))))
 494
 495 #+sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 496 (defoptimizer (log derive-type) ((x &optional y))
 497   (when (and (csubtypep (lvar-type x)
 498                         (specifier-type '(real 0.0)))
 499              (or (null y)
 500                  (csubtypep (lvar-type y)
 501                             (specifier-type '(real 0.0)))))
 502     (specifier-type 'float)))
 503 \f
 504 ;;;; irrational transforms
 505
 506 (defknown (%tan %sinh %asinh %atanh %log %logb %log10 %tan-quick)
 507           (double-float) double-float
 508   (movable foldable flushable))
 509
 510 (defknown (%sin %cos %tanh %sin-quick %cos-quick)
 511   (double-float) (double-float -1.0d0 1.0d0)
 512   (movable foldable flushable))
 513
 514 (defknown (%asin %atan)
 515   (double-float)
 516   (double-float #.(coerce (- (/ pi 2)) 'double-float)
 517                 #.(coerce (/ pi 2) 'double-float))
 518   (movable foldable flushable))
 519
 520 (defknown (%acos)
 521   (double-float) (double-float 0.0d0 #.(coerce pi 'double-float))
 522   (movable foldable flushable))
 523
 524 (defknown (%cosh)
 525   (double-float) (double-float 1.0d0)
 526   (movable foldable flushable))
 527
 528 (defknown (%acosh %exp %sqrt)
 529   (double-float) (double-float 0.0d0)
 530   (movable foldable flushable))
 531
 532 (defknown %expm1
 533   (double-float) (double-float -1d0)
 534   (movable foldable flushable))
 535
 536 (defknown (%hypot)
 537   (double-float double-float) (double-float 0d0)
 538   (movable foldable flushable))
 539
 540 (defknown (%pow)
 541   (double-float double-float) double-float
 542   (movable foldable flushable))
 543
 544 (defknown (%atan2)
 545   (double-float double-float)
 546   (double-float #.(coerce (- pi) 'double-float)
 547                 #.(coerce pi 'double-float))
 548   (movable foldable flushable))
 549
 550 (defknown (%scalb)
 551   (double-float double-float) double-float
 552   (movable foldable flushable))
 553
 554 (defknown (%scalbn)
 555   (double-float (signed-byte 32)) double-float
 556   (movable foldable flushable))
 557
 558 (defknown (%log1p)
 559   (double-float) double-float
 560   (movable foldable flushable))
 561
 562 (macrolet ((def (name prim rtype)
 563              `(progn
 564                (deftransform ,name ((x) (single-float) ,rtype)
 565                  `(coerce (,',prim (coerce x 'double-float)) 'single-float))
 566                (deftransform ,name ((x) (double-float) ,rtype)
 567                  `(,',prim x)))))
 568   (def exp %exp *)
 569   (def log %log float)
 570   (def sqrt %sqrt float)
 571   (def asin %asin float)
 572   (def acos %acos float)
 573   (def atan %atan *)
 574   (def sinh %sinh *)
 575   (def cosh %cosh *)
 576   (def tanh %tanh *)
 577   (def asinh %asinh *)
 578   (def acosh %acosh float)
 579   (def atanh %atanh float))
 580
 581 ;;; The argument range is limited on the x86 FP trig. functions. A
 582 ;;; post-test can detect a failure (and load a suitable result), but
 583 ;;; this test is avoided if possible.
 584 (macrolet ((def (name prim prim-quick)
 585              (declare (ignorable prim-quick))
 586              `(progn
 587                 (deftransform ,name ((x) (single-float) *)
 588                   #!+x86 (cond ((csubtypep (lvar-type x)
 589                                            (specifier-type '(single-float
 590                                                              (#.(- (expt 2f0 63)))
 591                                                              (#.(expt 2f0 63)))))
 592                                 `(coerce (,',prim-quick (coerce x 'double-float))
 593                                   'single-float))
 594                                (t
 595                                 (compiler-notify
 596                                  "unable to avoid inline argument range check~@
 597                                   because the argument range (~S) was not within 2^63"
 598                                  (type-specifier (lvar-type x)))
 599                                 `(coerce (,',prim (coerce x 'double-float)) 'single-float)))
 600                   #!-x86 `(coerce (,',prim (coerce x 'double-float)) 'single-float))
 601                (deftransform ,name ((x) (double-float) *)
 602                  #!+x86 (cond ((csubtypep (lvar-type x)
 603                                           (specifier-type '(double-float
 604                                                             (#.(- (expt 2d0 63)))
 605                                                             (#.(expt 2d0 63)))))
 606                                `(,',prim-quick x))
 607                               (t
 608                                (compiler-notify
 609                                 "unable to avoid inline argument range check~@
 610                                  because the argument range (~S) was not within 2^63"
 611                                 (type-specifier (lvar-type x)))
 612                                `(,',prim x)))
 613                  #!-x86 `(,',prim x)))))
 614   (def sin %sin %sin-quick)
 615   (def cos %cos %cos-quick)
 616   (def tan %tan %tan-quick))
 617
 618 (deftransform atan ((x y) (single-float single-float) *)
 619   `(coerce (%atan2 (coerce x 'double-float) (coerce y 'double-float))
 620     'single-float))
 621 (deftransform atan ((x y) (double-float double-float) *)
 622   `(%atan2 x y))
 623
 624 (deftransform expt ((x y) ((single-float 0f0) single-float) *)
 625   `(coerce (%pow (coerce x 'double-float) (coerce y 'double-float))
 626     'single-float))
 627 (deftransform expt ((x y) ((double-float 0d0) double-float) *)
 628   `(%pow x y))
 629 (deftransform expt ((x y) ((single-float 0f0) (signed-byte 32)) *)
 630   `(coerce (%pow (coerce x 'double-float) (coerce y 'double-float))
 631     'single-float))
 632 (deftransform expt ((x y) ((double-float 0d0) (signed-byte 32)) *)
 633   `(%pow x (coerce y 'double-float)))
 634
 635 ;;; ANSI says log with base zero returns zero.
 636 (deftransform log ((x y) (float float) float)
 637   '(if (zerop y) y (/ (log x) (log y))))
 638 \f
 639 ;;; Handle some simple transformations.
 640
 641 (deftransform abs ((x) ((complex double-float)) double-float)
 642   '(%hypot (realpart x) (imagpart x)))
 643
 644 (deftransform abs ((x) ((complex single-float)) single-float)
 645   '(coerce (%hypot (coerce (realpart x) 'double-float)
 646                    (coerce (imagpart x) 'double-float))
 647           'single-float))
 648
 649 (deftransform phase ((x) ((complex double-float)) double-float)
 650   '(%atan2 (imagpart x) (realpart x)))
 651
 652 (deftransform phase ((x) ((complex single-float)) single-float)
 653   '(coerce (%atan2 (coerce (imagpart x) 'double-float)
 654                    (coerce (realpart x) 'double-float))
 655           'single-float))
 656
 657 (deftransform phase ((x) ((float)) float)
 658   '(if (minusp (float-sign x))
 659        (float pi x)
 660        (float 0 x)))
 661
 662 ;;; The number is of type REAL.
 663 (defun numeric-type-real-p (type)
 664   (and (numeric-type-p type)
 665        (eq (numeric-type-complexp type) :real)))
 666
 667 ;;; Coerce a numeric type bound to the given type while handling
 668 ;;; exclusive bounds.
 669 (defun coerce-numeric-bound (bound type)
 670   (when bound
 671     (if (consp bound)
 672         (list (coerce (car bound) type))
 673         (coerce bound type))))
 674
 675 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 676 (progn
 677
 678 ;;;; optimizers for elementary functions
 679 ;;;;
 680 ;;;; These optimizers compute the output range of the elementary
 681 ;;;; function, based on the domain of the input.
 682
 683 ;;; Generate a specifier for a complex type specialized to the same
 684 ;;; type as the argument.
 685 (defun complex-float-type (arg)
 686   (declare (type numeric-type arg))
 687   (let* ((format (case (numeric-type-class arg)
 688                    ((integer rational) 'single-float)
 689                    (t (numeric-type-format arg))))
 690          (float-type (or format 'float)))
 691     (specifier-type `(complex ,float-type))))
 692
 693 ;;; Compute a specifier like '(OR FLOAT (COMPLEX FLOAT)), except float
 694 ;;; should be the right kind of float. Allow bounds for the float
 695 ;;; part too.
 696 (defun float-or-complex-float-type (arg &optional lo hi)
 697   (declare (type numeric-type arg))
 698   (let* ((format (case (numeric-type-class arg)
 699                    ((integer rational) 'single-float)
 700                    (t (numeric-type-format arg))))
 701          (float-type (or format 'float))
 702          (lo (coerce-numeric-bound lo float-type))
 703          (hi (coerce-numeric-bound hi float-type)))
 704     (specifier-type `(or (,float-type ,(or lo '*) ,(or hi '*))
 705                          (complex ,float-type)))))
 706
 707 ) ; PROGN
 708
 709 (eval-when (:compile-toplevel :execute)
 710   ;; So the problem with this hack is that it's actually broken.  If
 711   ;; the host does not have long floats, then setting *R-D-F-F* to
 712   ;; LONG-FLOAT doesn't actually buy us anything.  FIXME.
 713   (setf *read-default-float-format*
 714         #!+long-float 'long-float #!-long-float 'double-float))
 715 ;;; Test whether the numeric-type ARG is within the domain specified by
 716 ;;; DOMAIN-LOW and DOMAIN-HIGH, consider negative and positive zero to
 717 ;;; be distinct.
 718 #-sb-xc-host  ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 719 (defun domain-subtypep (arg domain-low domain-high)
 720   (declare (type numeric-type arg)
 721            (type (or real null) domain-low domain-high))
 722   (let* ((arg-lo (numeric-type-low arg))
 723          (arg-lo-val (type-bound-number arg-lo))
 724          (arg-hi (numeric-type-high arg))
 725          (arg-hi-val (type-bound-number arg-hi)))
 726     ;; Check that the ARG bounds are correctly canonicalized.
 727     (when (and arg-lo (floatp arg-lo-val) (zerop arg-lo-val) (consp arg-lo)
 728                (minusp (float-sign arg-lo-val)))
 729       (compiler-notify "float zero bound ~S not correctly canonicalized?" arg-lo)
 730       (setq arg-lo 0e0 arg-lo-val arg-lo))
 731     (when (and arg-hi (zerop arg-hi-val) (floatp arg-hi-val) (consp arg-hi)
 732                (plusp (float-sign arg-hi-val)))
 733       (compiler-notify "float zero bound ~S not correctly canonicalized?" arg-hi)
 734       (setq arg-hi (ecase *read-default-float-format*
 735                      (double-float (load-time-value (make-unportable-float :double-float-negative-zero)))
 736                      #!+long-float
 737                      (long-float (load-time-value (make-unportable-float :long-float-negative-zero))))
 738             arg-hi-val arg-hi))
 739     (flet ((fp-neg-zero-p (f)           ; Is F -0.0?
 740              (and (floatp f) (zerop f) (minusp (float-sign f))))
 741            (fp-pos-zero-p (f)           ; Is F +0.0?
 742              (and (floatp f) (zerop f) (plusp (float-sign f)))))
 743       (and (or (null domain-low)
 744                (and arg-lo (>= arg-lo-val domain-low)
 745                     (not (and (fp-pos-zero-p domain-low)
 746                               (fp-neg-zero-p arg-lo)))))
 747            (or (null domain-high)
 748                (and arg-hi (<= arg-hi-val domain-high)
 749                     (not (and (fp-neg-zero-p domain-high)
 750                               (fp-pos-zero-p arg-hi)))))))))
 751 (eval-when (:compile-toplevel :execute)
 752   (setf *read-default-float-format* 'single-float))
 753
 754 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 755 (progn
 756
 757 ;;; Handle monotonic functions of a single variable whose domain is
 758 ;;; possibly part of the real line. ARG is the variable, FUN is the
 759 ;;; function, and DOMAIN is a specifier that gives the (real) domain
 760 ;;; of the function. If ARG is a subset of the DOMAIN, we compute the
 761 ;;; bounds directly. Otherwise, we compute the bounds for the
 762 ;;; intersection between ARG and DOMAIN, and then append a complex
 763 ;;; result, which occurs for the parts of ARG not in the DOMAIN.
 764 ;;;
 765 ;;; Negative and positive zero are considered distinct within
 766 ;;; DOMAIN-LOW and DOMAIN-HIGH.
 767 ;;;
 768 ;;; DEFAULT-LOW and DEFAULT-HIGH are the lower and upper bounds if we
 769 ;;; can't compute the bounds using FUN.
 770 (defun elfun-derive-type-simple (arg fun domain-low domain-high
 771                                      default-low default-high
 772                                      &optional (increasingp t))
 773   (declare (type (or null real) domain-low domain-high))
 774   (etypecase arg
 775     (numeric-type
 776      (cond ((eq (numeric-type-complexp arg) :complex)
 777             (complex-float-type arg))
 778            ((numeric-type-real-p arg)
 779             ;; The argument is real, so let's find the intersection
 780             ;; between the argument and the domain of the function.
 781             ;; We compute the bounds on the intersection, and for
 782             ;; everything else, we return a complex number of the
 783             ;; appropriate type.
 784             (multiple-value-bind (intersection difference)
 785                 (interval-intersection/difference (numeric-type->interval arg)
 786                                                   (make-interval
 787                                                    :low domain-low
 788                                                    :high domain-high))
 789               (cond
 790                 (intersection
 791                  ;; Process the intersection.
 792                  (let* ((low (interval-low intersection))
 793                         (high (interval-high intersection))
 794                         (res-lo (or (bound-func fun (if increasingp low high))
 795                                     default-low))
 796                         (res-hi (or (bound-func fun (if increasingp high low))
 797                                     default-high))
 798                         (format (case (numeric-type-class arg)
 799                                   ((integer rational) 'single-float)
 800                                   (t (numeric-type-format arg))))
 801                         (bound-type (or format 'float))
 802                         (result-type
 803                          (make-numeric-type
 804                           :class 'float
 805                           :format format
 806                           :low (coerce-numeric-bound res-lo bound-type)
 807                           :high (coerce-numeric-bound res-hi bound-type))))
 808                    ;; If the ARG is a subset of the domain, we don't
 809                    ;; have to worry about the difference, because that
 810                    ;; can't occur.
 811                    (if (or (null difference)
 812                            ;; Check whether the arg is within the domain.
 813                            (domain-subtypep arg domain-low domain-high))
 814                        result-type
 815                        (list result-type
 816                              (specifier-type `(complex ,bound-type))))))
 817                 (t
 818                  ;; No intersection so the result must be purely complex.
 819                  (complex-float-type arg)))))
 820            (t
 821             (float-or-complex-float-type arg default-low default-high))))))
 822
 823 (macrolet
 824     ((frob (name domain-low domain-high def-low-bnd def-high-bnd
 825                  &key (increasingp t))
 826        (let ((num (gensym)))
 827          `(defoptimizer (,name derive-type) ((,num))
 828            (one-arg-derive-type
 829             ,num
 830             (lambda (arg)
 831               (elfun-derive-type-simple arg #',name
 832                                         ,domain-low ,domain-high
 833                                         ,def-low-bnd ,def-high-bnd
 834                                         ,increasingp))
 835             #',name)))))
 836   ;; These functions are easy because they are defined for the whole
 837   ;; real line.
 838   (frob exp nil nil 0 nil)
 839   (frob sinh nil nil nil nil)
 840   (frob tanh nil nil -1 1)
 841   (frob asinh nil nil nil nil)
 842
 843   ;; These functions are only defined for part of the real line. The
 844   ;; condition selects the desired part of the line.
 845   (frob asin -1d0 1d0 (- (/ pi 2)) (/ pi 2))
 846   ;; Acos is monotonic decreasing, so we need to swap the function
 847   ;; values at the lower and upper bounds of the input domain.
 848   (frob acos -1d0 1d0 0 pi :increasingp nil)
 849   (frob acosh 1d0 nil nil nil)
 850   (frob atanh -1d0 1d0 -1 1)
 851   ;; Kahan says that (sqrt -0.0) is -0.0, so use a specifier that
 852   ;; includes -0.0.
 853   (frob sqrt (load-time-value (make-unportable-float :double-float-negative-zero)) nil 0 nil))
 854
 855 ;;; Compute bounds for (expt x y). This should be easy since (expt x
 856 ;;; y) = (exp (* y (log x))). However, computations done this way
 857 ;;; have too much roundoff. Thus we have to do it the hard way.
 858 (defun safe-expt (x y)
 859   (handler-case
 860       (when (< (abs y) 10000)
 861         (expt x y))
 862     (error ()
 863       nil)))
 864
 865 ;;; Handle the case when x >= 1.
 866 (defun interval-expt-> (x y)
 867   (case (sb!c::interval-range-info y 0d0)
 868     (+
 869      ;; Y is positive and log X >= 0. The range of exp(y * log(x)) is
 870      ;; obviously non-negative. We just have to be careful for
 871      ;; infinite bounds (given by nil).
 872      (let ((lo (safe-expt (type-bound-number (sb!c::interval-low x))
 873                           (type-bound-number (sb!c::interval-low y))))
 874            (hi (safe-expt (type-bound-number (sb!c::interval-high x))
 875                           (type-bound-number (sb!c::interval-high y)))))
 876        (list (sb!c::make-interval :low (or lo 1) :high hi))))
 877     (-
 878      ;; Y is negative and log x >= 0. The range of exp(y * log(x)) is
 879      ;; obviously [0, 1]. However, underflow (nil) means 0 is the
 880      ;; result.
 881      (let ((lo (safe-expt (type-bound-number (sb!c::interval-high x))
 882                           (type-bound-number (sb!c::interval-low y))))
 883            (hi (safe-expt (type-bound-number (sb!c::interval-low x))
 884                           (type-bound-number (sb!c::interval-high y)))))
 885        (list (sb!c::make-interval :low (or lo 0) :high (or hi 1)))))
 886     (t
 887      ;; Split the interval in half.
 888      (destructuring-bind (y- y+)
 889          (sb!c::interval-split 0 y t)
 890        (list (interval-expt-> x y-)
 891              (interval-expt-> x y+))))))
 892
 893 ;;; Handle the case when x <= 1
 894 (defun interval-expt-< (x y)
 895   (case (sb!c::interval-range-info x 0d0)
 896     (+
 897      ;; The case of 0 <= x <= 1 is easy
 898      (case (sb!c::interval-range-info y)
 899        (+
 900         ;; Y is positive and log X <= 0. The range of exp(y * log(x)) is
 901         ;; obviously [0, 1]. We just have to be careful for infinite bounds
 902         ;; (given by nil).
 903         (let ((lo (safe-expt (type-bound-number (sb!c::interval-low x))
 904                              (type-bound-number (sb!c::interval-high y))))
 905               (hi (safe-expt (type-bound-number (sb!c::interval-high x))
 906                              (type-bound-number (sb!c::interval-low y)))))
 907           (list (sb!c::make-interval :low (or lo 0) :high (or hi 1)))))
 908        (-
 909         ;; Y is negative and log x <= 0. The range of exp(y * log(x)) is
 910         ;; obviously [1, inf].
 911         (let ((hi (safe-expt (type-bound-number (sb!c::interval-low x))
 912                              (type-bound-number (sb!c::interval-low y))))
 913               (lo (safe-expt (type-bound-number (sb!c::interval-high x))
 914                              (type-bound-number (sb!c::interval-high y)))))
 915           (list (sb!c::make-interval :low (or lo 1) :high hi))))
 916        (t
 917         ;; Split the interval in half
 918         (destructuring-bind (y- y+)
 919             (sb!c::interval-split 0 y t)
 920           (list (interval-expt-< x y-)
 921                 (interval-expt-< x y+))))))
 922     (-
 923      ;; The case where x <= 0. Y MUST be an INTEGER for this to work!
 924      ;; The calling function must insure this! For now we'll just
 925      ;; return the appropriate unbounded float type.
 926      (list (sb!c::make-interval :low nil :high nil)))
 927     (t
 928      (destructuring-bind (neg pos)
 929          (interval-split 0 x t t)
 930        (list (interval-expt-< neg y)
 931              (interval-expt-< pos y))))))
 932
 933 ;;; Compute bounds for (expt x y).
 934 (defun interval-expt (x y)
 935   (case (interval-range-info x 1)
 936     (+
 937      ;; X >= 1
 938          (interval-expt-> x y))
 939     (-
 940      ;; X <= 1
 941      (interval-expt-< x y))
 942     (t
 943      (destructuring-bind (left right)
 944          (interval-split 1 x t t)
 945        (list (interval-expt left y)
 946              (interval-expt right y))))))
 947
 948 (defun fixup-interval-expt (bnd x-int y-int x-type y-type)
 949   (declare (ignore x-int))
 950   ;; Figure out what the return type should be, given the argument
 951   ;; types and bounds and the result type and bounds.
 952   (cond ((csubtypep x-type (specifier-type 'integer))
 953          ;; an integer to some power
 954          (case (numeric-type-class y-type)
 955            (integer
 956             ;; Positive integer to an integer power is either an
 957             ;; integer or a rational.
 958             (let ((lo (or (interval-low bnd) '*))
 959                   (hi (or (interval-high bnd) '*)))
 960               (if (and (interval-low y-int)
 961                        (>= (type-bound-number (interval-low y-int)) 0))
 962                   (specifier-type `(integer ,lo ,hi))
 963                   (specifier-type `(rational ,lo ,hi)))))
 964            (rational
 965             ;; Positive integer to rational power is either a rational
 966             ;; or a single-float.
 967             (let* ((lo (interval-low bnd))
 968                    (hi (interval-high bnd))
 969                    (int-lo (if lo
 970                                (floor (type-bound-number lo))
 971                                '*))
 972                    (int-hi (if hi
 973                                (ceiling (type-bound-number hi))
 974                                '*))
 975                    (f-lo (or (bound-func #'float lo)
 976                              '*))
 977                    (f-hi (or (bound-func #'float hi)
 978                              '*)))
 979               (specifier-type `(or (rational ,int-lo ,int-hi)
 980                                 (single-float ,f-lo, f-hi)))))
 981            (float
 982             ;; A positive integer to a float power is a float.
 983             (modified-numeric-type y-type
 984                                    :low (interval-low bnd)
 985                                    :high (interval-high bnd)))
 986            (t
 987             ;; A positive integer to a number is a number (for now).
 988             (specifier-type 'number))))
 989         ((csubtypep x-type (specifier-type 'rational))
 990          ;; a rational to some power
 991          (case (numeric-type-class y-type)
 992            (integer
 993             ;; A positive rational to an integer power is always a rational.
 994             (specifier-type `(rational ,(or (interval-low bnd) '*)
 995                                        ,(or (interval-high bnd) '*))))
 996            (rational
 997             ;; A positive rational to rational power is either a rational
 998             ;; or a single-float.
 999             (let* ((lo (interval-low bnd))
1000                    (hi (interval-high bnd))
1001                    (int-lo (if lo
1002                                (floor (type-bound-number lo))
1003                                '*))
1004                    (int-hi (if hi
1005                                (ceiling (type-bound-number hi))
1006                                '*))
1007                    (f-lo (or (bound-func #'float lo)
1008                              '*))
1009                    (f-hi (or (bound-func #'float hi)
1010                              '*)))
1011               (specifier-type `(or (rational ,int-lo ,int-hi)
1012                                 (single-float ,f-lo, f-hi)))))
1013            (float
1014             ;; A positive rational to a float power is a float.
1015             (modified-numeric-type y-type
1016                                    :low (interval-low bnd)
1017                                    :high (interval-high bnd)))
1018            (t
1019             ;; A positive rational to a number is a number (for now).
1020             (specifier-type 'number))))
1021         ((csubtypep x-type (specifier-type 'float))
1022          ;; a float to some power
1023          (case (numeric-type-class y-type)
1024            ((or integer rational)
1025             ;; A positive float to an integer or rational power is
1026             ;; always a float.
1027             (make-numeric-type
1028              :class 'float
1029              :format (numeric-type-format x-type)
1030              :low (interval-low bnd)
1031              :high (interval-high bnd)))
1032            (float
1033             ;; A positive float to a float power is a float of the
1034             ;; higher type.
1035             (make-numeric-type
1036              :class 'float
1037              :format (float-format-max (numeric-type-format x-type)
1038                                        (numeric-type-format y-type))
1039              :low (interval-low bnd)
1040              :high (interval-high bnd)))
1041            (t
1042             ;; A positive float to a number is a number (for now)
1043             (specifier-type 'number))))
1044         (t
1045          ;; A number to some power is a number.
1046          (specifier-type 'number))))
1047
1048 (defun merged-interval-expt (x y)
1049   (let* ((x-int (numeric-type->interval x))
1050          (y-int (numeric-type->interval y)))
1051     (mapcar (lambda (type)
1052               (fixup-interval-expt type x-int y-int x y))
1053             (flatten-list (interval-expt x-int y-int)))))
1054
1055 (defun expt-derive-type-aux (x y same-arg)
1056   (declare (ignore same-arg))
1057   (cond ((or (not (numeric-type-real-p x))
1058              (not (numeric-type-real-p y)))
1059          ;; Use numeric contagion if either is not real.
1060          (numeric-contagion x y))
1061         ((csubtypep y (specifier-type 'integer))
1062          ;; A real raised to an integer power is well-defined.
1063          (merged-interval-expt x y))
1064         ;; A real raised to a non-integral power can be a float or a
1065         ;; complex number.
1066         ((or (csubtypep x (specifier-type '(rational 0)))
1067              (csubtypep x (specifier-type '(float (0d0)))))
1068          ;; But a positive real to any power is well-defined.
1069          (merged-interval-expt x y))
1070         ((and (csubtypep x (specifier-type 'rational))
1071               (csubtypep y (specifier-type 'rational)))
1072          ;; A rational to the power of a rational could be a rational
1073          ;; or a possibly-complex single float
1074          (specifier-type '(or rational single-float (complex single-float))))
1075         (t
1076          ;; a real to some power. The result could be a real or a
1077          ;; complex.
1078          (float-or-complex-float-type (numeric-contagion x y)))))
1079
1080 (defoptimizer (expt derive-type) ((x y))
1081   (two-arg-derive-type x y #'expt-derive-type-aux #'expt))
1082
1083 ;;; Note we must assume that a type including 0.0 may also include
1084 ;;; -0.0 and thus the result may be complex -infinity + i*pi.
1085 (defun log-derive-type-aux-1 (x)
1086   (elfun-derive-type-simple x #'log 0d0 nil nil nil))
1087
1088 (defun log-derive-type-aux-2 (x y same-arg)
1089   (let ((log-x (log-derive-type-aux-1 x))
1090         (log-y (log-derive-type-aux-1 y))
1091         (accumulated-list nil))
1092     ;; LOG-X or LOG-Y might be union types. We need to run through
1093     ;; the union types ourselves because /-DERIVE-TYPE-AUX doesn't.
1094     (dolist (x-type (prepare-arg-for-derive-type log-x))
1095       (dolist (y-type (prepare-arg-for-derive-type log-y))
1096         (push (/-derive-type-aux x-type y-type same-arg) accumulated-list)))
1097     (apply #'type-union (flatten-list accumulated-list))))
1098
1099 (defoptimizer (log derive-type) ((x &optional y))
1100   (if y
1101       (two-arg-derive-type x y #'log-derive-type-aux-2 #'log)
1102       (one-arg-derive-type x #'log-derive-type-aux-1 #'log)))
1103
1104 (defun atan-derive-type-aux-1 (y)
1105   (elfun-derive-type-simple y #'atan nil nil (- (/ pi 2)) (/ pi 2)))
1106
1107 (defun atan-derive-type-aux-2 (y x same-arg)
1108   (declare (ignore same-arg))
1109   ;; The hard case with two args. We just return the max bounds.
1110   (let ((result-type (numeric-contagion y x)))
1111     (cond ((and (numeric-type-real-p x)
1112                 (numeric-type-real-p y))
1113            (let* (;; FIXME: This expression for FORMAT seems to
1114                   ;; appear multiple times, and should be factored out.
1115                   (format (case (numeric-type-class result-type)
1116                             ((integer rational) 'single-float)
1117                             (t (numeric-type-format result-type))))
1118                   (bound-format (or format 'float)))
1119              (make-numeric-type :class 'float
1120                                 :format format
1121                                 :complexp :real
1122                                 :low (coerce (- pi) bound-format)
1123                                 :high (coerce pi bound-format))))
1124           (t
1125            ;; The result is a float or a complex number
1126            (float-or-complex-float-type result-type)))))
1127
1128 (defoptimizer (atan derive-type) ((y &optional x))
1129   (if x
1130       (two-arg-derive-type y x #'atan-derive-type-aux-2 #'atan)
1131       (one-arg-derive-type y #'atan-derive-type-aux-1 #'atan)))
1132
1133 (defun cosh-derive-type-aux (x)
1134   ;; We note that cosh x = cosh |x| for all real x.
1135   (elfun-derive-type-simple
1136    (if (numeric-type-real-p x)
1137        (abs-derive-type-aux x)
1138        x)
1139    #'cosh nil nil 0 nil))
1140
1141 (defoptimizer (cosh derive-type) ((num))
1142   (one-arg-derive-type num #'cosh-derive-type-aux #'cosh))
1143
1144 (defun phase-derive-type-aux (arg)
1145   (let* ((format (case (numeric-type-class arg)
1146                    ((integer rational) 'single-float)
1147                    (t (numeric-type-format arg))))
1148          (bound-type (or format 'float)))
1149     (cond ((numeric-type-real-p arg)
1150            (case (interval-range-info (numeric-type->interval arg) 0.0)
1151              (+
1152               ;; The number is positive, so the phase is 0.
1153               (make-numeric-type :class 'float
1154                                  :format format
1155                                  :complexp :real
1156                                  :low (coerce 0 bound-type)
1157                                  :high (coerce 0 bound-type)))
1158              (-
1159               ;; The number is always negative, so the phase is pi.
1160               (make-numeric-type :class 'float
1161                                  :format format
1162                                  :complexp :real
1163                                  :low (coerce pi bound-type)
1164                                  :high (coerce pi bound-type)))
1165              (t
1166               ;; We can't tell. The result is 0 or pi. Use a union
1167               ;; type for this.
1168               (list
1169                (make-numeric-type :class 'float
1170                                   :format format
1171                                   :complexp :real
1172                                   :low (coerce 0 bound-type)
1173                                   :high (coerce 0 bound-type))
1174                (make-numeric-type :class 'float
1175                                   :format format
1176                                   :complexp :real
1177                                   :low (coerce pi bound-type)
1178                                   :high (coerce pi bound-type))))))
1179           (t
1180            ;; We have a complex number. The answer is the range -pi
1181            ;; to pi. (-pi is included because we have -0.)
1182            (make-numeric-type :class 'float
1183                               :format format
1184                               :complexp :real
1185                               :low (coerce (- pi) bound-type)
1186                               :high (coerce pi bound-type))))))
1187
1188 (defoptimizer (phase derive-type) ((num))
1189   (one-arg-derive-type num #'phase-derive-type-aux #'phase))
1190
1191 ) ; PROGN
1192
1193 (deftransform realpart ((x) ((complex rational)) *)
1194   '(sb!kernel:%realpart x))
1195 (deftransform imagpart ((x) ((complex rational)) *)
1196   '(sb!kernel:%imagpart x))
1197
1198 ;;; Make REALPART and IMAGPART return the appropriate types. This
1199 ;;; should help a lot in optimized code.
1200 (defun realpart-derive-type-aux (type)
1201   (let ((class (numeric-type-class type))
1202         (format (numeric-type-format type)))
1203     (cond ((numeric-type-real-p type)
1204            ;; The realpart of a real has the same type and range as
1205            ;; the input.
1206            (make-numeric-type :class class
1207                               :format format
1208                               :complexp :real
1209                               :low (numeric-type-low type)
1210                               :high (numeric-type-high type)))
1211           (t
1212            ;; We have a complex number. The result has the same type
1213            ;; as the real part, except that it's real, not complex,
1214            ;; obviously.
1215            (make-numeric-type :class class
1216                               :format format
1217                               :complexp :real
1218                               :low (numeric-type-low type)
1219                               :high (numeric-type-high type))))))
1220 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1221 (defoptimizer (realpart derive-type) ((num))
1222   (one-arg-derive-type num #'realpart-derive-type-aux #'realpart))
1223 (defun imagpart-derive-type-aux (type)
1224   (let ((class (numeric-type-class type))
1225         (format (numeric-type-format type)))
1226     (cond ((numeric-type-real-p type)
1227            ;; The imagpart of a real has the same type as the input,
1228            ;; except that it's zero.
1229            (let ((bound-format (or format class 'real)))
1230              (make-numeric-type :class class
1231                                 :format format
1232                                 :complexp :real
1233                                 :low (coerce 0 bound-format)
1234                                 :high (coerce 0 bound-format))))
1235           (t
1236            ;; We have a complex number. The result has the same type as
1237            ;; the imaginary part, except that it's real, not complex,
1238            ;; obviously.
1239            (make-numeric-type :class class
1240                               :format format
1241                               :complexp :real
1242                               :low (numeric-type-low type)
1243                               :high (numeric-type-high type))))))
1244 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1245 (defoptimizer (imagpart derive-type) ((num))
1246   (one-arg-derive-type num #'imagpart-derive-type-aux #'imagpart))
1247
1248 (defun complex-derive-type-aux-1 (re-type)
1249   (if (numeric-type-p re-type)
1250       (make-numeric-type :class (numeric-type-class re-type)
1251                          :format (numeric-type-format re-type)
1252                          :complexp (if (csubtypep re-type
1253                                                   (specifier-type 'rational))
1254                                        :real
1255                                        :complex)
1256                          :low (numeric-type-low re-type)
1257                          :high (numeric-type-high re-type))
1258       (specifier-type 'complex)))
1259
1260 (defun complex-derive-type-aux-2 (re-type im-type same-arg)
1261   (declare (ignore same-arg))
1262   (if (and (numeric-type-p re-type)
1263            (numeric-type-p im-type))
1264       ;; Need to check to make sure numeric-contagion returns the
1265       ;; right type for what we want here.
1266
1267       ;; Also, what about rational canonicalization, like (complex 5 0)
1268       ;; is 5?  So, if the result must be complex, we make it so.
1269       ;; If the result might be complex, which happens only if the
1270       ;; arguments are rational, we make it a union type of (or
1271       ;; rational (complex rational)).
1272       (let* ((element-type (numeric-contagion re-type im-type))
1273              (rat-result-p (csubtypep element-type
1274                                       (specifier-type 'rational))))
1275         (if rat-result-p
1276             (type-union element-type
1277                         (specifier-type
1278                          `(complex ,(numeric-type-class element-type))))
1279             (make-numeric-type :class (numeric-type-class element-type)
1280                                :format (numeric-type-format element-type)
1281                                :complexp (if rat-result-p
1282                                              :real
1283                                              :complex))))
1284       (specifier-type 'complex)))
1285
1286 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1287 (defoptimizer (complex derive-type) ((re &optional im))
1288   (if im
1289       (two-arg-derive-type re im #'complex-derive-type-aux-2 #'complex)
1290       (one-arg-derive-type re #'complex-derive-type-aux-1 #'complex)))
1291
1292 ;;; Define some transforms for complex operations. We do this in lieu
1293 ;;; of complex operation VOPs.
1294 (macrolet ((frob (type)
1295              `(progn
1296                 (deftransform complex ((r) (,type))
1297                   '(complex r ,(coerce 0 type)))
1298                 (deftransform complex ((r i) (,type (and real (not ,type))))
1299                   '(complex r (truly-the ,type (coerce i ',type))))
1300                 (deftransform complex ((r i) ((and real (not ,type)) ,type))
1301                   '(complex (truly-the ,type (coerce r ',type)) i))
1302                ;; negation
1303                 #!-complex-float-vops
1304                (deftransform %negate ((z) ((complex ,type)) *)
1305                  '(complex (%negate (realpart z)) (%negate (imagpart z))))
1306                ;; complex addition and subtraction
1307                #!-complex-float-vops
1308                (deftransform + ((w z) ((complex ,type) (complex ,type)) *)
1309                  '(complex (+ (realpart w) (realpart z))
1310                            (+ (imagpart w) (imagpart z))))
1311                #!-complex-float-vops
1312                (deftransform - ((w z) ((complex ,type) (complex ,type)) *)
1313                  '(complex (- (realpart w) (realpart z))
1314                            (- (imagpart w) (imagpart z))))
1315                ;; Add and subtract a complex and a real.
1316                #!-complex-float-vops
1317                (deftransform + ((w z) ((complex ,type) real) *)
1318                  `(complex (+ (realpart w) z)
1319                            (+ (imagpart w) ,(coerce 0 ',type))))
1320                #!-complex-float-vops
1321                (deftransform + ((z w) (real (complex ,type)) *)
1322                  `(complex (+ (realpart w) z)
1323                            (+ (imagpart w) ,(coerce 0 ',type))))
1324                ;; Add and subtract a real and a complex number.
1325                #!-complex-float-vops
1326                (deftransform - ((w z) ((complex ,type) real) *)
1327                  `(complex (- (realpart w) z)
1328                            (- (imagpart w) ,(coerce 0 ',type))))
1329                #!-complex-float-vops
1330                (deftransform - ((z w) (real (complex ,type)) *)
1331                  `(complex (- z (realpart w))
1332                            (- ,(coerce 0 ',type) (imagpart w))))
1333                ;; Multiply and divide two complex numbers.
1334                #!-complex-float-vops
1335                (deftransform * ((x y) ((complex ,type) (complex ,type)) *)
1336                  '(let* ((rx (realpart x))
1337                          (ix (imagpart x))
1338                          (ry (realpart y))
1339                          (iy (imagpart y)))
1340                     (complex (- (* rx ry) (* ix iy))
1341                              (+ (* rx iy) (* ix ry)))))
1342                (deftransform / ((x y) ((complex ,type) (complex ,type)) *)
1343                  #!-complex-float-vops
1344                  '(let* ((rx (realpart x))
1345                          (ix (imagpart x))
1346                          (ry (realpart y))
1347                          (iy (imagpart y)))
1348                     (if (> (abs ry) (abs iy))
1349                         (let* ((r (/ iy ry))
1350                                (dn (+ ry (* r iy))))
1351                           (complex (/ (+ rx (* ix r)) dn)
1352                                    (/ (- ix (* rx r)) dn)))
1353                         (let* ((r (/ ry iy))
1354                                (dn (+ iy (* r ry))))
1355                           (complex (/ (+ (* rx r) ix) dn)
1356                                    (/ (- (* ix r) rx) dn)))))
1357                  #!+complex-float-vops
1358                  `(let* ((cs (conjugate (sb!vm::swap-complex x)))
1359                          (ry (realpart y))
1360                          (iy (imagpart y)))
1361                     (if (> (abs ry) (abs iy))
1362                         (let* ((r (/ iy ry))
1363                                (dn (+ ry (* r iy))))
1364                           (/ (+ x (* cs r)) dn))
1365                         (let* ((r (/ ry iy))
1366                                (dn (+ iy (* r ry))))
1367                           (/ (+ (* x r) cs) dn)))))
1368                ;; Multiply a complex by a real or vice versa.
1369                #!-complex-float-vops
1370                (deftransform * ((w z) ((complex ,type) real) *)
1371                  '(complex (* (realpart w) z) (* (imagpart w) z)))
1372                #!-complex-float-vops
1373                (deftransform * ((z w) (real (complex ,type)) *)
1374                  '(complex (* (realpart w) z) (* (imagpart w) z)))
1375                ;; Divide a complex by a real or vice versa.
1376                #!-complex-float-vops
1377                (deftransform / ((w z) ((complex ,type) real) *)
1378                  '(complex (/ (realpart w) z) (/ (imagpart w) z)))
1379                (deftransform / ((x y) (,type (complex ,type)) *)
1380                  #!-complex-float-vops
1381                  '(let* ((ry (realpart y))
1382                          (iy (imagpart y)))
1383                     (if (> (abs ry) (abs iy))
1384                         (let* ((r (/ iy ry))
1385                                (dn (+ ry (* r iy))))
1386                           (complex (/ x dn)
1387                                    (/ (- (* x r)) dn)))
1388                         (let* ((r (/ ry iy))
1389                                (dn (+ iy (* r ry))))
1390                           (complex (/ (* x r) dn)
1391                                    (/ (- x) dn)))))
1392                  #!+complex-float-vops
1393                  '(let* ((ry (realpart y))
1394                          (iy (imagpart y)))
1395                    (if (> (abs ry) (abs iy))
1396                        (let* ((r (/ iy ry))
1397                               (dn (+ ry (* r iy))))
1398                          (/ (complex x (- (* x r))) dn))
1399                        (let* ((r (/ ry iy))
1400                               (dn (+ iy (* r ry))))
1401                          (/ (complex (* x r) (- x)) dn)))))
1402                ;; conjugate of complex number
1403                #!-complex-float-vops
1404                (deftransform conjugate ((z) ((complex ,type)) *)
1405                  '(complex (realpart z) (- (imagpart z))))
1406                ;; CIS
1407                (deftransform cis ((z) ((,type)) *)
1408                  '(complex (cos z) (sin z)))
1409                ;; comparison
1410                #!-complex-float-vops
1411                (deftransform = ((w z) ((complex ,type) (complex ,type)) *)
1412                  '(and (= (realpart w) (realpart z))
1413                        (= (imagpart w) (imagpart z))))
1414                #!-complex-float-vops
1415                (deftransform = ((w z) ((complex ,type) real) *)
1416                  '(and (= (realpart w) z) (zerop (imagpart w))))
1417                #!-complex-float-vops
1418                (deftransform = ((w z) (real (complex ,type)) *)
1419                  '(and (= (realpart z) w) (zerop (imagpart z)))))))
1420
1421   (frob single-float)
1422   (frob double-float))
1423
1424 ;;; Here are simple optimizers for SIN, COS, and TAN. They do not
1425 ;;; produce a minimal range for the result; the result is the widest
1426 ;;; possible answer. This gets around the problem of doing range
1427 ;;; reduction correctly but still provides useful results when the
1428 ;;; inputs are union types.
1429 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1430 (progn
1431 (defun trig-derive-type-aux (arg domain fun
1432                                  &optional def-lo def-hi (increasingp t))
1433   (etypecase arg
1434     (numeric-type
1435      (cond ((eq (numeric-type-complexp arg) :complex)
1436             (make-numeric-type :class (numeric-type-class arg)
1437                                :format (numeric-type-format arg)
1438                                :complexp :complex))
1439            ((numeric-type-real-p arg)
1440             (let* ((format (case (numeric-type-class arg)
1441                              ((integer rational) 'single-float)
1442                              (t (numeric-type-format arg))))
1443                    (bound-type (or format 'float)))
1444               ;; If the argument is a subset of the "principal" domain
1445               ;; of the function, we can compute the bounds because
1446               ;; the function is monotonic. We can't do this in
1447               ;; general for these periodic functions because we can't
1448               ;; (and don't want to) do the argument reduction in
1449               ;; exactly the same way as the functions themselves do
1450               ;; it.
1451               (if (csubtypep arg domain)
1452                   (let ((res-lo (bound-func fun (numeric-type-low arg)))
1453                         (res-hi (bound-func fun (numeric-type-high arg))))
1454                     (unless increasingp
1455                       (rotatef res-lo res-hi))
1456                     (make-numeric-type
1457                      :class 'float
1458                      :format format
1459                      :low (coerce-numeric-bound res-lo bound-type)
1460                      :high (coerce-numeric-bound res-hi bound-type)))
1461                   (make-numeric-type
1462                    :class 'float
1463                    :format format
1464                    :low (and def-lo (coerce def-lo bound-type))
1465                    :high (and def-hi (coerce def-hi bound-type))))))
1466            (t
1467             (float-or-complex-float-type arg def-lo def-hi))))))
1468
1469 (defoptimizer (sin derive-type) ((num))
1470   (one-arg-derive-type
1471    num
1472    (lambda (arg)
1473      ;; Derive the bounds if the arg is in [-pi/2, pi/2].
1474      (trig-derive-type-aux
1475       arg
1476       (specifier-type `(float ,(- (/ pi 2)) ,(/ pi 2)))
1477       #'sin
1478       -1 1))
1479    #'sin))
1480
1481 (defoptimizer (cos derive-type) ((num))
1482   (one-arg-derive-type
1483    num
1484    (lambda (arg)
1485      ;; Derive the bounds if the arg is in [0, pi].
1486      (trig-derive-type-aux arg
1487                            (specifier-type `(float 0d0 ,pi))
1488                            #'cos
1489                            -1 1
1490                            nil))
1491    #'cos))
1492
1493 (defoptimizer (tan derive-type) ((num))
1494   (one-arg-derive-type
1495    num
1496    (lambda (arg)
1497      ;; Derive the bounds if the arg is in [-pi/2, pi/2].
1498      (trig-derive-type-aux arg
1499                            (specifier-type `(float ,(- (/ pi 2)) ,(/ pi 2)))
1500                            #'tan
1501                            nil nil))
1502    #'tan))
1503
1504 (defoptimizer (conjugate derive-type) ((num))
1505   (one-arg-derive-type num
1506     (lambda (arg)
1507       (flet ((most-negative-bound (l h)
1508                (and l h
1509                     (if (< (type-bound-number l) (- (type-bound-number h)))
1510                         l
1511                         (set-bound (- (type-bound-number h)) (consp h)))))
1512              (most-positive-bound (l h)
1513                (and l h
1514                     (if (> (type-bound-number h) (- (type-bound-number l)))
1515                         h
1516                         (set-bound (- (type-bound-number l)) (consp l))))))
1517         (if (numeric-type-real-p arg)
1518             (lvar-type num)
1519             (let ((low (numeric-type-low arg))
1520                   (high (numeric-type-high arg)))
1521               (let ((new-low (most-negative-bound low high))
1522                     (new-high (most-positive-bound low high)))
1523               (modified-numeric-type arg :low new-low :high new-high))))))
1524     #'conjugate))
1525
1526 (defoptimizer (cis derive-type) ((num))
1527   (one-arg-derive-type num
1528     (lambda (arg)
1529       (sb!c::specifier-type
1530        `(complex ,(or (numeric-type-format arg) 'float))))
1531     #'cis))
1532
1533 ) ; PROGN
1534 \f
1535 ;;;; TRUNCATE, FLOOR, CEILING, and ROUND
1536
1537 (macrolet ((define-frobs (fun ufun)
1538              `(progn
1539                 (defknown ,ufun (real) integer (movable foldable flushable))
1540                 (deftransform ,fun ((x &optional by)
1541                                     (* &optional
1542                                        (constant-arg (member 1))))
1543                   '(let ((res (,ufun x)))
1544                      (values res (- x res)))))))
1545   (define-frobs truncate %unary-truncate)
1546   (define-frobs round %unary-round))
1547
1548 (deftransform %unary-truncate ((x) (single-float))
1549   `(%unary-truncate/single-float x))
1550 (deftransform %unary-truncate ((x) (double-float))
1551   `(%unary-truncate/double-float x))
1552
1553 ;;; Convert (TRUNCATE x y) to the obvious implementation.
1554 ;;;
1555 ;;; ...plus hair: Insert explicit coercions to appropriate float types: Python
1556 ;;; is reluctant it generate explicit integer->float coercions due to
1557 ;;; precision issues (see SAFE-SINGLE-COERCION-P &co), but this is not an
1558 ;;; issue here as there is no DERIVE-TYPE optimizer on specialized versions of
1559 ;;; %UNARY-TRUNCATE, so the derived type of TRUNCATE remains the best we can
1560 ;;; do here -- which is fine. Also take care not to add unnecassary division
1561 ;;; or multiplication by 1, since we are not able to always eliminate them,
1562 ;;; depending on FLOAT-ACCURACY. Finally, leave out the secondary value when
1563 ;;; we know it is unused: COERCE is not flushable.
1564 (macrolet ((def (type other-float-arg-types)
1565              (let ((unary (symbolicate "%UNARY-TRUNCATE/" type))
1566                    (coerce (symbolicate "%" type)))
1567                `(deftransform truncate ((x &optional y)
1568                                         (,type
1569                                          &optional (or ,type ,@other-float-arg-types integer))
1570                                         * :result result)
1571                   (let* ((result-type (and result
1572                                            (lvar-derived-type result)))
1573                          (compute-all (and (values-type-p result-type)
1574                                            (not (type-single-value-p result-type)))))
1575                     (if (or (not y)
1576                             (and (constant-lvar-p y) (= 1 (lvar-value y))))
1577                         (if compute-all
1578                             `(let ((res (,',unary x)))
1579                                (values res (- x (,',coerce res))))
1580                             `(let ((res (,',unary x)))
1581                                ;; Dummy secondary value!
1582                                (values res x)))
1583                         (if compute-all
1584                             `(let* ((f (,',coerce y))
1585                                     (res (,',unary (/ x f))))
1586                                (values res (- x (* f (,',coerce res)))))
1587                             `(let* ((f (,',coerce y))
1588                                     (res (,',unary (/ x f))))
1589                                ;; Dummy secondary value!
1590                                (values res x)))))))))
1591   (def single-float ())
1592   (def double-float (single-float)))
1593
1594 (deftransform floor ((number &optional divisor)
1595                      (float &optional (or integer float)))
1596   (let ((defaulted-divisor (if divisor 'divisor 1)))
1597     `(multiple-value-bind (tru rem) (truncate number ,defaulted-divisor)
1598        (if (and (not (zerop rem))
1599                 (if (minusp ,defaulted-divisor)
1600                     (plusp number)
1601                     (minusp number)))
1602            (values (1- tru) (+ rem ,defaulted-divisor))
1603            (values tru rem)))))
1604
1605 (deftransform ceiling ((number &optional divisor)
1606                        (float &optional (or integer float)))
1607   (let ((defaulted-divisor (if divisor 'divisor 1)))
1608     `(multiple-value-bind (tru rem) (truncate number ,defaulted-divisor)
1609        (if (and (not (zerop rem))
1610                 (if (minusp ,defaulted-divisor)
1611                     (minusp number)
1612                     (plusp number)))
1613            (values (1+ tru) (- rem ,defaulted-divisor))
1614            (values tru rem)))))
1615
1616 (defknown %unary-ftruncate (real) float (movable foldable flushable))
1617 (defknown %unary-ftruncate/single (single-float) single-float
1618   (movable foldable flushable))
1619 (defknown %unary-ftruncate/double (double-float) double-float
1620   (movable foldable flushable))
1621
1622 (defun %unary-ftruncate/single (x)
1623   (declare (type single-float x))
1624   (declare (optimize speed (safety 0)))
1625   (let* ((bits (single-float-bits x))
1626          (exp (ldb sb!vm:single-float-exponent-byte bits))
1627          (biased (the single-float-exponent
1628                    (- exp sb!vm:single-float-bias))))
1629     (declare (type (signed-byte 32) bits))
1630     (cond
1631       ((= exp sb!vm:single-float-normal-exponent-max) x)
1632       ((<= biased 0) (* x 0f0))
1633       ((>= biased (float-digits x)) x)
1634       (t
1635        (let ((frac-bits (- (float-digits x) biased)))
1636          (setf bits (logandc2 bits (- (ash 1 frac-bits) 1)))
1637          (make-single-float bits))))))
1638
1639 (defun %unary-ftruncate/double (x)
1640   (declare (type double-float x))
1641   (declare (optimize speed (safety 0)))
1642   (let* ((high (double-float-high-bits x))
1643          (low (double-float-low-bits x))
1644          (exp (ldb sb!vm:double-float-exponent-byte high))
1645          (biased (the double-float-exponent
1646                    (- exp sb!vm:double-float-bias))))
1647     (declare (type (signed-byte 32) high)
1648              (type (unsigned-byte 32) low))
1649     (cond
1650       ((= exp sb!vm:double-float-normal-exponent-max) x)
1651       ((<= biased 0) (* x 0d0))
1652       ((>= biased (float-digits x)) x)
1653       (t
1654        (let ((frac-bits (- (float-digits x) biased)))
1655          (cond ((< frac-bits 32)
1656                 (setf low (logandc2 low (- (ash 1 frac-bits) 1))))
1657                (t
1658                 (setf low 0)
1659                 (setf high (logandc2 high (- (ash 1 (- frac-bits 32)) 1)))))
1660          (make-double-float high low))))))
1661
1662 (macrolet
1663     ((def (float-type fun)
1664          `(deftransform %unary-ftruncate ((x) (,float-type))
1665             '(,fun x))))
1666   (def single-float %unary-ftruncate/single)
1667   (def double-float %unary-ftruncate/double))