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                                      nil))
 331                      (hi (bound-func (lambda (x)
 332                                        (if (< ,most-positive x )
 333                                            ,most-positive
 334                                            (coerce x ',type)))
 335                                      (numeric-type-high num)
 336                                      nil)))
 337                 (specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
 338
 339             (defoptimizer (,fun derive-type) ((num))
 340               (handler-case
 341                   (one-arg-derive-type num #',aux-name #',fun)
 342                 (type-error ()
 343                   nil)))))))
 344   (frob %single-float single-float
 345         most-negative-single-float most-positive-single-float)
 346   (frob %double-float double-float
 347         most-negative-double-float most-positive-double-float))
 348 ) ; PROGN
 349 \f
 350 ;;;; float contagion
 351
 352 (defun safe-ctype-for-single-coercion-p (x)
 353   ;; See comment in SAFE-SINGLE-COERCION-P -- this deals with the same
 354   ;; problem, but in the context of evaluated and compiled (+ <int> <single>)
 355   ;; giving different result if we fail to check for this.
 356   (or (not (csubtypep x (specifier-type 'integer)))
 357       #!+x86
 358       (csubtypep x (specifier-type `(integer ,most-negative-exactly-single-float-fixnum
 359                                              ,most-positive-exactly-single-float-fixnum)))
 360       #!-x86
 361       (csubtypep x (specifier-type 'fixnum))))
 362
 363 ;;; Do some stuff to recognize when the loser is doing mixed float and
 364 ;;; rational arithmetic, or different float types, and fix it up. If
 365 ;;; we don't, he won't even get so much as an efficiency note.
 366 (deftransform float-contagion-arg1 ((x y) * * :defun-only t :node node)
 367   (if (or (not (types-equal-or-intersect (lvar-type y) (specifier-type 'single-float)))
 368           (safe-ctype-for-single-coercion-p (lvar-type x)))
 369       `(,(lvar-fun-name (basic-combination-fun node))
 370          (float x y) y)
 371       (give-up-ir1-transform)))
 372 (deftransform float-contagion-arg2 ((x y) * * :defun-only t :node node)
 373   (if (or (not (types-equal-or-intersect (lvar-type x) (specifier-type 'single-float)))
 374           (safe-ctype-for-single-coercion-p (lvar-type y)))
 375       `(,(lvar-fun-name (basic-combination-fun node))
 376          x (float y x))
 377       (give-up-ir1-transform)))
 378
 379 (dolist (x '(+ * / -))
 380   (%deftransform x '(function (rational float) *) #'float-contagion-arg1)
 381   (%deftransform x '(function (float rational) *) #'float-contagion-arg2))
 382
 383 (dolist (x '(= < > + * / -))
 384   (%deftransform x '(function (single-float double-float) *)
 385                  #'float-contagion-arg1)
 386   (%deftransform x '(function (double-float single-float) *)
 387                  #'float-contagion-arg2))
 388
 389 (macrolet ((def (type &rest args)
 390              `(deftransform * ((x y) (,type (constant-arg (member ,@args))) *
 391                                ;; Beware the SNaN!
 392                                :policy (zerop float-accuracy))
 393                 "optimize multiplication by one"
 394                 (let ((y (lvar-value y)))
 395                   (if (minusp y)
 396                       '(%negate x)
 397                       'x)))))
 398   (def single-float 1.0 -1.0)
 399   (def double-float 1.0d0 -1.0d0))
 400
 401 ;;; Return the reciprocal of X if it can be represented exactly, NIL otherwise.
 402 (defun maybe-exact-reciprocal (x)
 403   (unless (zerop x)
 404     (handler-case
 405         (multiple-value-bind (significand exponent sign)
 406             (integer-decode-float x)
 407           ;; only powers of 2 can be inverted exactly
 408           (unless (zerop (logand significand (1- significand)))
 409             (return-from maybe-exact-reciprocal nil))
 410           (let ((expected   (/ sign significand (expt 2 exponent)))
 411                 (reciprocal (/ x)))
 412             (multiple-value-bind (significand exponent sign)
 413                 (integer-decode-float reciprocal)
 414               ;; Denorms can't be inverted safely.
 415               (and (eql expected (* sign significand (expt 2 exponent)))
 416                    reciprocal))))
 417       (error () (return-from maybe-exact-reciprocal nil)))))
 418
 419 ;;; Replace constant division by multiplication with exact reciprocal,
 420 ;;; if one exists.
 421 (macrolet ((def (type)
 422              `(deftransform / ((x y) (,type (constant-arg ,type)) *
 423                                :node node)
 424                 "convert to multiplication by reciprocal"
 425                 (let ((n (lvar-value y)))
 426                   (if (policy node (zerop float-accuracy))
 427                       `(* x ,(/ n))
 428                       (let ((r (maybe-exact-reciprocal n)))
 429                         (if r
 430                             `(* x ,r)
 431                             (give-up-ir1-transform
 432                              "~S does not have an exact reciprocal"
 433                              n))))))))
 434   (def single-float)
 435   (def double-float))
 436
 437 ;;; Optimize addition and subtraction of zero
 438 (macrolet ((def (op type &rest args)
 439              `(deftransform ,op ((x y) (,type (constant-arg (member ,@args))) *
 440                                  ;; Beware the SNaN!
 441                                  :policy (zerop float-accuracy))
 442                 'x)))
 443   ;; No signed zeros, thanks.
 444   (def + single-float 0 0.0)
 445   (def - single-float 0 0.0)
 446   (def + double-float 0 0.0 0.0d0)
 447   (def - double-float 0 0.0 0.0d0))
 448
 449 ;;; On most platforms (+ x x) is faster than (* x 2)
 450 (macrolet ((def (type &rest args)
 451              `(deftransform * ((x y) (,type (constant-arg (member ,@args))))
 452                 '(+ x x))))
 453   (def single-float 2 2.0)
 454   (def double-float 2 2.0 2.0d0))
 455
 456 ;;; Prevent ZEROP, PLUSP, and MINUSP from losing horribly. We can't in
 457 ;;; general float rational args to comparison, since Common Lisp
 458 ;;; semantics says we are supposed to compare as rationals, but we can
 459 ;;; do it for any rational that has a precise representation as a
 460 ;;; float (such as 0).
 461 (macrolet ((frob (op)
 462              `(deftransform ,op ((x y) (float rational) *)
 463                 "open-code FLOAT to RATIONAL comparison"
 464                 (unless (constant-lvar-p y)
 465                   (give-up-ir1-transform
 466                    "The RATIONAL value isn't known at compile time."))
 467                 (let ((val (lvar-value y)))
 468                   (unless (eql (rational (float val)) val)
 469                     (give-up-ir1-transform
 470                      "~S doesn't have a precise float representation."
 471                      val)))
 472                 `(,',op x (float y x)))))
 473   (frob <)
 474   (frob >)
 475   (frob =))
 476 \f
 477 ;;;; irrational derive-type methods
 478
 479 ;;; Derive the result to be float for argument types in the
 480 ;;; appropriate domain.
 481 #+sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 482 (dolist (stuff '((asin (real -1.0 1.0))
 483                  (acos (real -1.0 1.0))
 484                  (acosh (real 1.0))
 485                  (atanh (real -1.0 1.0))
 486                  (sqrt (real 0.0))))
 487   (destructuring-bind (name type) stuff
 488     (let ((type (specifier-type type)))
 489       (setf (fun-info-derive-type (fun-info-or-lose name))
 490             (lambda (call)
 491               (declare (type combination call))
 492               (when (csubtypep (lvar-type
 493                                 (first (combination-args call)))
 494                                type)
 495                 (specifier-type 'float)))))))
 496
 497 #+sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 498 (defoptimizer (log derive-type) ((x &optional y))
 499   (when (and (csubtypep (lvar-type x)
 500                         (specifier-type '(real 0.0)))
 501              (or (null y)
 502                  (csubtypep (lvar-type y)
 503                             (specifier-type '(real 0.0)))))
 504     (specifier-type 'float)))
 505 \f
 506 ;;;; irrational transforms
 507
 508 (defknown (%tan %sinh %asinh %atanh %log %logb %log10 %tan-quick)
 509           (double-float) double-float
 510   (movable foldable flushable))
 511
 512 (defknown (%sin %cos %tanh %sin-quick %cos-quick)
 513   (double-float) (double-float -1.0d0 1.0d0)
 514   (movable foldable flushable))
 515
 516 (defknown (%asin %atan)
 517   (double-float)
 518   (double-float #.(coerce (- (/ pi 2)) 'double-float)
 519                 #.(coerce (/ pi 2) 'double-float))
 520   (movable foldable flushable))
 521
 522 (defknown (%acos)
 523   (double-float) (double-float 0.0d0 #.(coerce pi 'double-float))
 524   (movable foldable flushable))
 525
 526 (defknown (%cosh)
 527   (double-float) (double-float 1.0d0)
 528   (movable foldable flushable))
 529
 530 (defknown (%acosh %exp %sqrt)
 531   (double-float) (double-float 0.0d0)
 532   (movable foldable flushable))
 533
 534 (defknown %expm1
 535   (double-float) (double-float -1d0)
 536   (movable foldable flushable))
 537
 538 (defknown (%hypot)
 539   (double-float double-float) (double-float 0d0)
 540   (movable foldable flushable))
 541
 542 (defknown (%pow)
 543   (double-float double-float) double-float
 544   (movable foldable flushable))
 545
 546 (defknown (%atan2)
 547   (double-float double-float)
 548   (double-float #.(coerce (- pi) 'double-float)
 549                 #.(coerce pi 'double-float))
 550   (movable foldable flushable))
 551
 552 (defknown (%scalb)
 553   (double-float double-float) double-float
 554   (movable foldable flushable))
 555
 556 (defknown (%scalbn)
 557   (double-float (signed-byte 32)) double-float
 558   (movable foldable flushable))
 559
 560 (defknown (%log1p)
 561   (double-float) double-float
 562   (movable foldable flushable))
 563
 564 (macrolet ((def (name prim rtype)
 565              `(progn
 566                (deftransform ,name ((x) (single-float) ,rtype)
 567                  `(coerce (,',prim (coerce x 'double-float)) 'single-float))
 568                (deftransform ,name ((x) (double-float) ,rtype)
 569                  `(,',prim x)))))
 570   (def exp %exp *)
 571   (def log %log float)
 572   (def sqrt %sqrt float)
 573   (def asin %asin float)
 574   (def acos %acos float)
 575   (def atan %atan *)
 576   (def sinh %sinh *)
 577   (def cosh %cosh *)
 578   (def tanh %tanh *)
 579   (def asinh %asinh *)
 580   (def acosh %acosh float)
 581   (def atanh %atanh float))
 582
 583 ;;; The argument range is limited on the x86 FP trig. functions. A
 584 ;;; post-test can detect a failure (and load a suitable result), but
 585 ;;; this test is avoided if possible.
 586 (macrolet ((def (name prim prim-quick)
 587              (declare (ignorable prim-quick))
 588              `(progn
 589                 (deftransform ,name ((x) (single-float) *)
 590                   #!+x86 (cond ((csubtypep (lvar-type x)
 591                                            (specifier-type '(single-float
 592                                                              (#.(- (expt 2f0 63)))
 593                                                              (#.(expt 2f0 63)))))
 594                                 `(coerce (,',prim-quick (coerce x 'double-float))
 595                                   'single-float))
 596                                (t
 597                                 (compiler-notify
 598                                  "unable to avoid inline argument range check~@
 599                                   because the argument range (~S) was not within 2^63"
 600                                  (type-specifier (lvar-type x)))
 601                                 `(coerce (,',prim (coerce x 'double-float)) 'single-float)))
 602                   #!-x86 `(coerce (,',prim (coerce x 'double-float)) 'single-float))
 603                (deftransform ,name ((x) (double-float) *)
 604                  #!+x86 (cond ((csubtypep (lvar-type x)
 605                                           (specifier-type '(double-float
 606                                                             (#.(- (expt 2d0 63)))
 607                                                             (#.(expt 2d0 63)))))
 608                                `(,',prim-quick x))
 609                               (t
 610                                (compiler-notify
 611                                 "unable to avoid inline argument range check~@
 612                                  because the argument range (~S) was not within 2^63"
 613                                 (type-specifier (lvar-type x)))
 614                                `(,',prim x)))
 615                  #!-x86 `(,',prim x)))))
 616   (def sin %sin %sin-quick)
 617   (def cos %cos %cos-quick)
 618   (def tan %tan %tan-quick))
 619
 620 (deftransform atan ((x y) (single-float single-float) *)
 621   `(coerce (%atan2 (coerce x 'double-float) (coerce y 'double-float))
 622     'single-float))
 623 (deftransform atan ((x y) (double-float double-float) *)
 624   `(%atan2 x y))
 625
 626 (deftransform expt ((x y) ((single-float 0f0) single-float) *)
 627   `(coerce (%pow (coerce x 'double-float) (coerce y 'double-float))
 628     'single-float))
 629 (deftransform expt ((x y) ((double-float 0d0) double-float) *)
 630   `(%pow x y))
 631 (deftransform expt ((x y) ((single-float 0f0) (signed-byte 32)) *)
 632   `(coerce (%pow (coerce x 'double-float) (coerce y 'double-float))
 633     'single-float))
 634 (deftransform expt ((x y) ((double-float 0d0) (signed-byte 32)) *)
 635   `(%pow x (coerce y 'double-float)))
 636
 637 ;;; ANSI says log with base zero returns zero.
 638 (deftransform log ((x y) (float float) float)
 639   '(if (zerop y) y (/ (log x) (log y))))
 640 \f
 641 ;;; Handle some simple transformations.
 642
 643 (deftransform abs ((x) ((complex double-float)) double-float)
 644   '(%hypot (realpart x) (imagpart x)))
 645
 646 (deftransform abs ((x) ((complex single-float)) single-float)
 647   '(coerce (%hypot (coerce (realpart x) 'double-float)
 648                    (coerce (imagpart x) 'double-float))
 649           'single-float))
 650
 651 (deftransform phase ((x) ((complex double-float)) double-float)
 652   '(%atan2 (imagpart x) (realpart x)))
 653
 654 (deftransform phase ((x) ((complex single-float)) single-float)
 655   '(coerce (%atan2 (coerce (imagpart x) 'double-float)
 656                    (coerce (realpart x) 'double-float))
 657           'single-float))
 658
 659 (deftransform phase ((x) ((float)) float)
 660   '(if (minusp (float-sign x))
 661        (float pi x)
 662        (float 0 x)))
 663
 664 ;;; The number is of type REAL.
 665 (defun numeric-type-real-p (type)
 666   (and (numeric-type-p type)
 667        (eq (numeric-type-complexp type) :real)))
 668
 669 ;;; Coerce a numeric type bound to the given type while handling
 670 ;;; exclusive bounds.
 671 (defun coerce-numeric-bound (bound type)
 672   (when bound
 673     (if (consp bound)
 674         (list (coerce (car bound) type))
 675         (coerce bound type))))
 676
 677 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 678 (progn
 679
 680 ;;;; optimizers for elementary functions
 681 ;;;;
 682 ;;;; These optimizers compute the output range of the elementary
 683 ;;;; function, based on the domain of the input.
 684
 685 ;;; Generate a specifier for a complex type specialized to the same
 686 ;;; type as the argument.
 687 (defun complex-float-type (arg)
 688   (declare (type numeric-type arg))
 689   (let* ((format (case (numeric-type-class arg)
 690                    ((integer rational) 'single-float)
 691                    (t (numeric-type-format arg))))
 692          (float-type (or format 'float)))
 693     (specifier-type `(complex ,float-type))))
 694
 695 ;;; Compute a specifier like '(OR FLOAT (COMPLEX FLOAT)), except float
 696 ;;; should be the right kind of float. Allow bounds for the float
 697 ;;; part too.
 698 (defun float-or-complex-float-type (arg &optional lo hi)
 699   (declare (type numeric-type arg))
 700   (let* ((format (case (numeric-type-class arg)
 701                    ((integer rational) 'single-float)
 702                    (t (numeric-type-format arg))))
 703          (float-type (or format 'float))
 704          (lo (coerce-numeric-bound lo float-type))
 705          (hi (coerce-numeric-bound hi float-type)))
 706     (specifier-type `(or (,float-type ,(or lo '*) ,(or hi '*))
 707                          (complex ,float-type)))))
 708
 709 ) ; PROGN
 710
 711 (eval-when (:compile-toplevel :execute)
 712   ;; So the problem with this hack is that it's actually broken.  If
 713   ;; the host does not have long floats, then setting *R-D-F-F* to
 714   ;; LONG-FLOAT doesn't actually buy us anything.  FIXME.
 715   (setf *read-default-float-format*
 716         #!+long-float 'long-float #!-long-float 'double-float))
 717 ;;; Test whether the numeric-type ARG is within the domain specified by
 718 ;;; DOMAIN-LOW and DOMAIN-HIGH, consider negative and positive zero to
 719 ;;; be distinct.
 720 #-sb-xc-host  ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 721 (defun domain-subtypep (arg domain-low domain-high)
 722   (declare (type numeric-type arg)
 723            (type (or real null) domain-low domain-high))
 724   (let* ((arg-lo (numeric-type-low arg))
 725          (arg-lo-val (type-bound-number arg-lo))
 726          (arg-hi (numeric-type-high arg))
 727          (arg-hi-val (type-bound-number arg-hi)))
 728     ;; Check that the ARG bounds are correctly canonicalized.
 729     (when (and arg-lo (floatp arg-lo-val) (zerop arg-lo-val) (consp arg-lo)
 730                (minusp (float-sign arg-lo-val)))
 731       (compiler-notify "float zero bound ~S not correctly canonicalized?" arg-lo)
 732       (setq arg-lo 0e0 arg-lo-val arg-lo))
 733     (when (and arg-hi (zerop arg-hi-val) (floatp arg-hi-val) (consp arg-hi)
 734                (plusp (float-sign arg-hi-val)))
 735       (compiler-notify "float zero bound ~S not correctly canonicalized?" arg-hi)
 736       (setq arg-hi (ecase *read-default-float-format*
 737                      (double-float (load-time-value (make-unportable-float :double-float-negative-zero)))
 738                      #!+long-float
 739                      (long-float (load-time-value (make-unportable-float :long-float-negative-zero))))
 740             arg-hi-val arg-hi))
 741     (flet ((fp-neg-zero-p (f)           ; Is F -0.0?
 742              (and (floatp f) (zerop f) (minusp (float-sign f))))
 743            (fp-pos-zero-p (f)           ; Is F +0.0?
 744              (and (floatp f) (zerop f) (plusp (float-sign f)))))
 745       (and (or (null domain-low)
 746                (and arg-lo (>= arg-lo-val domain-low)
 747                     (not (and (fp-pos-zero-p domain-low)
 748                               (fp-neg-zero-p arg-lo)))))
 749            (or (null domain-high)
 750                (and arg-hi (<= arg-hi-val domain-high)
 751                     (not (and (fp-neg-zero-p domain-high)
 752                               (fp-pos-zero-p arg-hi)))))))))
 753 (eval-when (:compile-toplevel :execute)
 754   (setf *read-default-float-format* 'single-float))
 755
 756 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
 757 (progn
 758
 759 ;;; Handle monotonic functions of a single variable whose domain is
 760 ;;; possibly part of the real line. ARG is the variable, FUN is the
 761 ;;; function, and DOMAIN is a specifier that gives the (real) domain
 762 ;;; of the function. If ARG is a subset of the DOMAIN, we compute the
 763 ;;; bounds directly. Otherwise, we compute the bounds for the
 764 ;;; intersection between ARG and DOMAIN, and then append a complex
 765 ;;; result, which occurs for the parts of ARG not in the DOMAIN.
 766 ;;;
 767 ;;; Negative and positive zero are considered distinct within
 768 ;;; DOMAIN-LOW and DOMAIN-HIGH.
 769 ;;;
 770 ;;; DEFAULT-LOW and DEFAULT-HIGH are the lower and upper bounds if we
 771 ;;; can't compute the bounds using FUN.
 772 (defun elfun-derive-type-simple (arg fun domain-low domain-high
 773                                      default-low default-high
 774                                      &optional (increasingp t))
 775   (declare (type (or null real) domain-low domain-high))
 776   (etypecase arg
 777     (numeric-type
 778      (cond ((eq (numeric-type-complexp arg) :complex)
 779             (complex-float-type arg))
 780            ((numeric-type-real-p arg)
 781             ;; The argument is real, so let's find the intersection
 782             ;; between the argument and the domain of the function.
 783             ;; We compute the bounds on the intersection, and for
 784             ;; everything else, we return a complex number of the
 785             ;; appropriate type.
 786             (multiple-value-bind (intersection difference)
 787                 (interval-intersection/difference (numeric-type->interval arg)
 788                                                   (make-interval
 789                                                    :low domain-low
 790                                                    :high domain-high))
 791               (cond
 792                 (intersection
 793                  ;; Process the intersection.
 794                  (let* ((low (interval-low intersection))
 795                         (high (interval-high intersection))
 796                         (res-lo (or (bound-func fun (if increasingp low high) nil)
 797                                     default-low))
 798                         (res-hi (or (bound-func fun (if increasingp high low) nil)
 799                                     default-high))
 800                         (format (case (numeric-type-class arg)
 801                                   ((integer rational) 'single-float)
 802                                   (t (numeric-type-format arg))))
 803                         (bound-type (or format 'float))
 804                         (result-type
 805                          (make-numeric-type
 806                           :class 'float
 807                           :format format
 808                           :low (coerce-numeric-bound res-lo bound-type)
 809                           :high (coerce-numeric-bound res-hi bound-type))))
 810                    ;; If the ARG is a subset of the domain, we don't
 811                    ;; have to worry about the difference, because that
 812                    ;; can't occur.
 813                    (if (or (null difference)
 814                            ;; Check whether the arg is within the domain.
 815                            (domain-subtypep arg domain-low domain-high))
 816                        result-type
 817                        (list result-type
 818                              (specifier-type `(complex ,bound-type))))))
 819                 (t
 820                  ;; No intersection so the result must be purely complex.
 821                  (complex-float-type arg)))))
 822            (t
 823             (float-or-complex-float-type arg default-low default-high))))))
 824
 825 (macrolet
 826     ((frob (name domain-low domain-high def-low-bnd def-high-bnd
 827                  &key (increasingp t))
 828        (let ((num (gensym)))
 829          `(defoptimizer (,name derive-type) ((,num))
 830            (one-arg-derive-type
 831             ,num
 832             (lambda (arg)
 833               (elfun-derive-type-simple arg #',name
 834                                         ,domain-low ,domain-high
 835                                         ,def-low-bnd ,def-high-bnd
 836                                         ,increasingp))
 837             #',name)))))
 838   ;; These functions are easy because they are defined for the whole
 839   ;; real line.
 840   (frob exp nil nil 0 nil)
 841   (frob sinh nil nil nil nil)
 842   (frob tanh nil nil -1 1)
 843   (frob asinh nil nil nil nil)
 844
 845   ;; These functions are only defined for part of the real line. The
 846   ;; condition selects the desired part of the line.
 847   (frob asin -1d0 1d0 (- (/ pi 2)) (/ pi 2))
 848   ;; Acos is monotonic decreasing, so we need to swap the function
 849   ;; values at the lower and upper bounds of the input domain.
 850   (frob acos -1d0 1d0 0 pi :increasingp nil)
 851   (frob acosh 1d0 nil nil nil)
 852   (frob atanh -1d0 1d0 -1 1)
 853   ;; Kahan says that (sqrt -0.0) is -0.0, so use a specifier that
 854   ;; includes -0.0.
 855   (frob sqrt (load-time-value (make-unportable-float :double-float-negative-zero)) nil 0 nil))
 856
 857 ;;; Compute bounds for (expt x y). This should be easy since (expt x
 858 ;;; y) = (exp (* y (log x))). However, computations done this way
 859 ;;; have too much roundoff. Thus we have to do it the hard way.
 860 (defun safe-expt (x y)
 861   (handler-case
 862       (when (< (abs y) 10000)
 863         (expt x y))
 864     (error ()
 865       nil)))
 866
 867 ;;; Handle the case when x >= 1.
 868 (defun interval-expt-> (x y)
 869   (case (sb!c::interval-range-info y 0d0)
 870     (+
 871      ;; Y is positive and log X >= 0. The range of exp(y * log(x)) is
 872      ;; obviously non-negative. We just have to be careful for
 873      ;; infinite bounds (given by nil).
 874      (let ((lo (safe-expt (type-bound-number (sb!c::interval-low x))
 875                           (type-bound-number (sb!c::interval-low y))))
 876            (hi (safe-expt (type-bound-number (sb!c::interval-high x))
 877                           (type-bound-number (sb!c::interval-high y)))))
 878        (list (sb!c::make-interval :low (or lo 1) :high hi))))
 879     (-
 880      ;; Y is negative and log x >= 0. The range of exp(y * log(x)) is
 881      ;; obviously [0, 1]. However, underflow (nil) means 0 is the
 882      ;; result.
 883      (let ((lo (safe-expt (type-bound-number (sb!c::interval-high x))
 884                           (type-bound-number (sb!c::interval-low y))))
 885            (hi (safe-expt (type-bound-number (sb!c::interval-low x))
 886                           (type-bound-number (sb!c::interval-high y)))))
 887        (list (sb!c::make-interval :low (or lo 0) :high (or hi 1)))))
 888     (t
 889      ;; Split the interval in half.
 890      (destructuring-bind (y- y+)
 891          (sb!c::interval-split 0 y t)
 892        (list (interval-expt-> x y-)
 893              (interval-expt-> x y+))))))
 894
 895 ;;; Handle the case when x <= 1
 896 (defun interval-expt-< (x y)
 897   (case (sb!c::interval-range-info x 0d0)
 898     (+
 899      ;; The case of 0 <= x <= 1 is easy
 900      (case (sb!c::interval-range-info y)
 901        (+
 902         ;; Y is positive and log X <= 0. The range of exp(y * log(x)) is
 903         ;; obviously [0, 1]. We just have to be careful for infinite bounds
 904         ;; (given by nil).
 905         (let ((lo (safe-expt (type-bound-number (sb!c::interval-low x))
 906                              (type-bound-number (sb!c::interval-high y))))
 907               (hi (safe-expt (type-bound-number (sb!c::interval-high x))
 908                              (type-bound-number (sb!c::interval-low y)))))
 909           (list (sb!c::make-interval :low (or lo 0) :high (or hi 1)))))
 910        (-
 911         ;; Y is negative and log x <= 0. The range of exp(y * log(x)) is
 912         ;; obviously [1, inf].
 913         (let ((hi (safe-expt (type-bound-number (sb!c::interval-low x))
 914                              (type-bound-number (sb!c::interval-low y))))
 915               (lo (safe-expt (type-bound-number (sb!c::interval-high x))
 916                              (type-bound-number (sb!c::interval-high y)))))
 917           (list (sb!c::make-interval :low (or lo 1) :high hi))))
 918        (t
 919         ;; Split the interval in half
 920         (destructuring-bind (y- y+)
 921             (sb!c::interval-split 0 y t)
 922           (list (interval-expt-< x y-)
 923                 (interval-expt-< x y+))))))
 924     (-
 925      ;; The case where x <= 0. Y MUST be an INTEGER for this to work!
 926      ;; The calling function must insure this! For now we'll just
 927      ;; return the appropriate unbounded float type.
 928      (list (sb!c::make-interval :low nil :high nil)))
 929     (t
 930      (destructuring-bind (neg pos)
 931          (interval-split 0 x t t)
 932        (list (interval-expt-< neg y)
 933              (interval-expt-< pos y))))))
 934
 935 ;;; Compute bounds for (expt x y).
 936 (defun interval-expt (x y)
 937   (case (interval-range-info x 1)
 938     (+
 939      ;; X >= 1
 940          (interval-expt-> x y))
 941     (-
 942      ;; X <= 1
 943      (interval-expt-< x y))
 944     (t
 945      (destructuring-bind (left right)
 946          (interval-split 1 x t t)
 947        (list (interval-expt left y)
 948              (interval-expt right y))))))
 949
 950 (defun fixup-interval-expt (bnd x-int y-int x-type y-type)
 951   (declare (ignore x-int))
 952   ;; Figure out what the return type should be, given the argument
 953   ;; types and bounds and the result type and bounds.
 954   (cond ((csubtypep x-type (specifier-type 'integer))
 955          ;; an integer to some power
 956          (case (numeric-type-class y-type)
 957            (integer
 958             ;; Positive integer to an integer power is either an
 959             ;; integer or a rational.
 960             (let ((lo (or (interval-low bnd) '*))
 961                   (hi (or (interval-high bnd) '*)))
 962               (if (and (interval-low y-int)
 963                        (>= (type-bound-number (interval-low y-int)) 0))
 964                   (specifier-type `(integer ,lo ,hi))
 965                   (specifier-type `(rational ,lo ,hi)))))
 966            (rational
 967             ;; Positive integer to rational power is either a rational
 968             ;; or a single-float.
 969             (let* ((lo (interval-low bnd))
 970                    (hi (interval-high bnd))
 971                    (int-lo (if lo
 972                                (floor (type-bound-number lo))
 973                                '*))
 974                    (int-hi (if hi
 975                                (ceiling (type-bound-number hi))
 976                                '*))
 977                    (f-lo (or (bound-func #'float lo nil)
 978                              '*))
 979                    (f-hi (or (bound-func #'float hi nil)
 980                              '*)))
 981               (specifier-type `(or (rational ,int-lo ,int-hi)
 982                                 (single-float ,f-lo, f-hi)))))
 983            (float
 984             ;; A positive integer to a float power is a float.
 985             (let ((format (numeric-type-format y-type)))
 986               (aver format)
 987               (modified-numeric-type
 988                y-type
 989                :low (coerce-numeric-bound (interval-low bnd) format)
 990                :high (coerce-numeric-bound (interval-high bnd) format))))
 991            (t
 992             ;; A positive integer to a number is a number (for now).
 993             (specifier-type 'number))))
 994         ((csubtypep x-type (specifier-type 'rational))
 995          ;; a rational to some power
 996          (case (numeric-type-class y-type)
 997            (integer
 998             ;; A positive rational to an integer power is always a rational.
 999             (specifier-type `(rational ,(or (interval-low bnd) '*)
1000                                        ,(or (interval-high bnd) '*))))
1001            (rational
1002             ;; A positive rational to rational power is either a rational
1003             ;; or a single-float.
1004             (let* ((lo (interval-low bnd))
1005                    (hi (interval-high bnd))
1006                    (int-lo (if lo
1007                                (floor (type-bound-number lo))
1008                                '*))
1009                    (int-hi (if hi
1010                                (ceiling (type-bound-number hi))
1011                                '*))
1012                    (f-lo (or (bound-func #'float lo nil)
1013                              '*))
1014                    (f-hi (or (bound-func #'float hi nil)
1015                              '*)))
1016               (specifier-type `(or (rational ,int-lo ,int-hi)
1017                                 (single-float ,f-lo, f-hi)))))
1018            (float
1019             ;; A positive rational to a float power is a float.
1020             (let ((format (numeric-type-format y-type)))
1021               (aver format)
1022               (modified-numeric-type
1023                y-type
1024                :low (coerce-numeric-bound (interval-low bnd) format)
1025                :high (coerce-numeric-bound (interval-high bnd) format))))
1026            (t
1027             ;; A positive rational to a number is a number (for now).
1028             (specifier-type 'number))))
1029         ((csubtypep x-type (specifier-type 'float))
1030          ;; a float to some power
1031          (case (numeric-type-class y-type)
1032            ((or integer rational)
1033             ;; A positive float to an integer or rational power is
1034             ;; always a float.
1035             (let ((format (numeric-type-format x-type)))
1036               (aver format)
1037               (make-numeric-type
1038                :class 'float
1039                :format format
1040                :low (coerce-numeric-bound (interval-low bnd) format)
1041                :high (coerce-numeric-bound (interval-high bnd) format))))
1042            (float
1043             ;; A positive float to a float power is a float of the
1044             ;; higher type.
1045             (let ((format (float-format-max (numeric-type-format x-type)
1046                                             (numeric-type-format y-type))))
1047               (aver format)
1048               (make-numeric-type
1049                :class 'float
1050                :format format
1051                :low (coerce-numeric-bound (interval-low bnd) format)
1052                :high (coerce-numeric-bound (interval-high bnd) format))))
1053            (t
1054             ;; A positive float to a number is a number (for now)
1055             (specifier-type 'number))))
1056         (t
1057          ;; A number to some power is a number.
1058          (specifier-type 'number))))
1059
1060 (defun merged-interval-expt (x y)
1061   (let* ((x-int (numeric-type->interval x))
1062          (y-int (numeric-type->interval y)))
1063     (mapcar (lambda (type)
1064               (fixup-interval-expt type x-int y-int x y))
1065             (flatten-list (interval-expt x-int y-int)))))
1066
1067 (defun expt-derive-type-aux (x y same-arg)
1068   (declare (ignore same-arg))
1069   (cond ((or (not (numeric-type-real-p x))
1070              (not (numeric-type-real-p y)))
1071          ;; Use numeric contagion if either is not real.
1072          (numeric-contagion x y))
1073         ((csubtypep y (specifier-type 'integer))
1074          ;; A real raised to an integer power is well-defined.
1075          (merged-interval-expt x y))
1076         ;; A real raised to a non-integral power can be a float or a
1077         ;; complex number.
1078         ((or (csubtypep x (specifier-type '(rational 0)))
1079              (csubtypep x (specifier-type '(float (0d0)))))
1080          ;; But a positive real to any power is well-defined.
1081          (merged-interval-expt x y))
1082         ((and (csubtypep x (specifier-type 'rational))
1083               (csubtypep y (specifier-type 'rational)))
1084          ;; A rational to the power of a rational could be a rational
1085          ;; or a possibly-complex single float
1086          (specifier-type '(or rational single-float (complex single-float))))
1087         (t
1088          ;; a real to some power. The result could be a real or a
1089          ;; complex.
1090          (float-or-complex-float-type (numeric-contagion x y)))))
1091
1092 (defoptimizer (expt derive-type) ((x y))
1093   (two-arg-derive-type x y #'expt-derive-type-aux #'expt))
1094
1095 ;;; Note we must assume that a type including 0.0 may also include
1096 ;;; -0.0 and thus the result may be complex -infinity + i*pi.
1097 (defun log-derive-type-aux-1 (x)
1098   (elfun-derive-type-simple x #'log 0d0 nil nil nil))
1099
1100 (defun log-derive-type-aux-2 (x y same-arg)
1101   (let ((log-x (log-derive-type-aux-1 x))
1102         (log-y (log-derive-type-aux-1 y))
1103         (accumulated-list nil))
1104     ;; LOG-X or LOG-Y might be union types. We need to run through
1105     ;; the union types ourselves because /-DERIVE-TYPE-AUX doesn't.
1106     (dolist (x-type (prepare-arg-for-derive-type log-x))
1107       (dolist (y-type (prepare-arg-for-derive-type log-y))
1108         (push (/-derive-type-aux x-type y-type same-arg) accumulated-list)))
1109     (apply #'type-union (flatten-list accumulated-list))))
1110
1111 (defoptimizer (log derive-type) ((x &optional y))
1112   (if y
1113       (two-arg-derive-type x y #'log-derive-type-aux-2 #'log)
1114       (one-arg-derive-type x #'log-derive-type-aux-1 #'log)))
1115
1116 (defun atan-derive-type-aux-1 (y)
1117   (elfun-derive-type-simple y #'atan nil nil (- (/ pi 2)) (/ pi 2)))
1118
1119 (defun atan-derive-type-aux-2 (y x same-arg)
1120   (declare (ignore same-arg))
1121   ;; The hard case with two args. We just return the max bounds.
1122   (let ((result-type (numeric-contagion y x)))
1123     (cond ((and (numeric-type-real-p x)
1124                 (numeric-type-real-p y))
1125            (let* (;; FIXME: This expression for FORMAT seems to
1126                   ;; appear multiple times, and should be factored out.
1127                   (format (case (numeric-type-class result-type)
1128                             ((integer rational) 'single-float)
1129                             (t (numeric-type-format result-type))))
1130                   (bound-format (or format 'float)))
1131              (make-numeric-type :class 'float
1132                                 :format format
1133                                 :complexp :real
1134                                 :low (coerce (- pi) bound-format)
1135                                 :high (coerce pi bound-format))))
1136           (t
1137            ;; The result is a float or a complex number
1138            (float-or-complex-float-type result-type)))))
1139
1140 (defoptimizer (atan derive-type) ((y &optional x))
1141   (if x
1142       (two-arg-derive-type y x #'atan-derive-type-aux-2 #'atan)
1143       (one-arg-derive-type y #'atan-derive-type-aux-1 #'atan)))
1144
1145 (defun cosh-derive-type-aux (x)
1146   ;; We note that cosh x = cosh |x| for all real x.
1147   (elfun-derive-type-simple
1148    (if (numeric-type-real-p x)
1149        (abs-derive-type-aux x)
1150        x)
1151    #'cosh nil nil 0 nil))
1152
1153 (defoptimizer (cosh derive-type) ((num))
1154   (one-arg-derive-type num #'cosh-derive-type-aux #'cosh))
1155
1156 (defun phase-derive-type-aux (arg)
1157   (let* ((format (case (numeric-type-class arg)
1158                    ((integer rational) 'single-float)
1159                    (t (numeric-type-format arg))))
1160          (bound-type (or format 'float)))
1161     (cond ((numeric-type-real-p arg)
1162            (case (interval-range-info (numeric-type->interval arg) 0.0)
1163              (+
1164               ;; The number is positive, so the phase is 0.
1165               (make-numeric-type :class 'float
1166                                  :format format
1167                                  :complexp :real
1168                                  :low (coerce 0 bound-type)
1169                                  :high (coerce 0 bound-type)))
1170              (-
1171               ;; The number is always negative, so the phase is pi.
1172               (make-numeric-type :class 'float
1173                                  :format format
1174                                  :complexp :real
1175                                  :low (coerce pi bound-type)
1176                                  :high (coerce pi bound-type)))
1177              (t
1178               ;; We can't tell. The result is 0 or pi. Use a union
1179               ;; type for this.
1180               (list
1181                (make-numeric-type :class 'float
1182                                   :format format
1183                                   :complexp :real
1184                                   :low (coerce 0 bound-type)
1185                                   :high (coerce 0 bound-type))
1186                (make-numeric-type :class 'float
1187                                   :format format
1188                                   :complexp :real
1189                                   :low (coerce pi bound-type)
1190                                   :high (coerce pi bound-type))))))
1191           (t
1192            ;; We have a complex number. The answer is the range -pi
1193            ;; to pi. (-pi is included because we have -0.)
1194            (make-numeric-type :class 'float
1195                               :format format
1196                               :complexp :real
1197                               :low (coerce (- pi) bound-type)
1198                               :high (coerce pi bound-type))))))
1199
1200 (defoptimizer (phase derive-type) ((num))
1201   (one-arg-derive-type num #'phase-derive-type-aux #'phase))
1202
1203 ) ; PROGN
1204
1205 (deftransform realpart ((x) ((complex rational)) *)
1206   '(sb!kernel:%realpart x))
1207 (deftransform imagpart ((x) ((complex rational)) *)
1208   '(sb!kernel:%imagpart x))
1209
1210 ;;; Make REALPART and IMAGPART return the appropriate types. This
1211 ;;; should help a lot in optimized code.
1212 (defun realpart-derive-type-aux (type)
1213   (let ((class (numeric-type-class type))
1214         (format (numeric-type-format type)))
1215     (cond ((numeric-type-real-p type)
1216            ;; The realpart of a real has the same type and range as
1217            ;; the input.
1218            (make-numeric-type :class class
1219                               :format format
1220                               :complexp :real
1221                               :low (numeric-type-low type)
1222                               :high (numeric-type-high type)))
1223           (t
1224            ;; We have a complex number. The result has the same type
1225            ;; as the real part, except that it's real, not complex,
1226            ;; obviously.
1227            (make-numeric-type :class class
1228                               :format format
1229                               :complexp :real
1230                               :low (numeric-type-low type)
1231                               :high (numeric-type-high type))))))
1232 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1233 (defoptimizer (realpart derive-type) ((num))
1234   (one-arg-derive-type num #'realpart-derive-type-aux #'realpart))
1235 (defun imagpart-derive-type-aux (type)
1236   (let ((class (numeric-type-class type))
1237         (format (numeric-type-format type)))
1238     (cond ((numeric-type-real-p type)
1239            ;; The imagpart of a real has the same type as the input,
1240            ;; except that it's zero.
1241            (let ((bound-format (or format class 'real)))
1242              (make-numeric-type :class class
1243                                 :format format
1244                                 :complexp :real
1245                                 :low (coerce 0 bound-format)
1246                                 :high (coerce 0 bound-format))))
1247           (t
1248            ;; We have a complex number. The result has the same type as
1249            ;; the imaginary part, except that it's real, not complex,
1250            ;; obviously.
1251            (make-numeric-type :class class
1252                               :format format
1253                               :complexp :real
1254                               :low (numeric-type-low type)
1255                               :high (numeric-type-high type))))))
1256 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1257 (defoptimizer (imagpart derive-type) ((num))
1258   (one-arg-derive-type num #'imagpart-derive-type-aux #'imagpart))
1259
1260 (defun complex-derive-type-aux-1 (re-type)
1261   (if (numeric-type-p re-type)
1262       (make-numeric-type :class (numeric-type-class re-type)
1263                          :format (numeric-type-format re-type)
1264                          :complexp (if (csubtypep re-type
1265                                                   (specifier-type 'rational))
1266                                        :real
1267                                        :complex)
1268                          :low (numeric-type-low re-type)
1269                          :high (numeric-type-high re-type))
1270       (specifier-type 'complex)))
1271
1272 (defun complex-derive-type-aux-2 (re-type im-type same-arg)
1273   (declare (ignore same-arg))
1274   (if (and (numeric-type-p re-type)
1275            (numeric-type-p im-type))
1276       ;; Need to check to make sure numeric-contagion returns the
1277       ;; right type for what we want here.
1278
1279       ;; Also, what about rational canonicalization, like (complex 5 0)
1280       ;; is 5?  So, if the result must be complex, we make it so.
1281       ;; If the result might be complex, which happens only if the
1282       ;; arguments are rational, we make it a union type of (or
1283       ;; rational (complex rational)).
1284       (let* ((element-type (numeric-contagion re-type im-type))
1285              (rat-result-p (csubtypep element-type
1286                                       (specifier-type 'rational))))
1287         (if rat-result-p
1288             (type-union element-type
1289                         (specifier-type
1290                          `(complex ,(numeric-type-class element-type))))
1291             (make-numeric-type :class (numeric-type-class element-type)
1292                                :format (numeric-type-format element-type)
1293                                :complexp (if rat-result-p
1294                                              :real
1295                                              :complex))))
1296       (specifier-type 'complex)))
1297
1298 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1299 (defoptimizer (complex derive-type) ((re &optional im))
1300   (if im
1301       (two-arg-derive-type re im #'complex-derive-type-aux-2 #'complex)
1302       (one-arg-derive-type re #'complex-derive-type-aux-1 #'complex)))
1303
1304 ;;; Define some transforms for complex operations. We do this in lieu
1305 ;;; of complex operation VOPs.
1306 (macrolet ((frob (type)
1307              `(progn
1308                 (deftransform complex ((r) (,type))
1309                   '(complex r ,(coerce 0 type)))
1310                 (deftransform complex ((r i) (,type (and real (not ,type))))
1311                   '(complex r (truly-the ,type (coerce i ',type))))
1312                 (deftransform complex ((r i) ((and real (not ,type)) ,type))
1313                   '(complex (truly-the ,type (coerce r ',type)) i))
1314                ;; negation
1315                 #!-complex-float-vops
1316                (deftransform %negate ((z) ((complex ,type)) *)
1317                  '(complex (%negate (realpart z)) (%negate (imagpart z))))
1318                ;; complex addition and subtraction
1319                #!-complex-float-vops
1320                (deftransform + ((w z) ((complex ,type) (complex ,type)) *)
1321                  '(complex (+ (realpart w) (realpart z))
1322                            (+ (imagpart w) (imagpart z))))
1323                #!-complex-float-vops
1324                (deftransform - ((w z) ((complex ,type) (complex ,type)) *)
1325                  '(complex (- (realpart w) (realpart z))
1326                            (- (imagpart w) (imagpart z))))
1327                ;; Add and subtract a complex and a real.
1328                #!-complex-float-vops
1329                (deftransform + ((w z) ((complex ,type) real) *)
1330                  `(complex (+ (realpart w) z)
1331                            (+ (imagpart w) ,(coerce 0 ',type))))
1332                #!-complex-float-vops
1333                (deftransform + ((z w) (real (complex ,type)) *)
1334                  `(complex (+ (realpart w) z)
1335                            (+ (imagpart w) ,(coerce 0 ',type))))
1336                ;; Add and subtract a real and a complex number.
1337                #!-complex-float-vops
1338                (deftransform - ((w z) ((complex ,type) real) *)
1339                  `(complex (- (realpart w) z)
1340                            (- (imagpart w) ,(coerce 0 ',type))))
1341                #!-complex-float-vops
1342                (deftransform - ((z w) (real (complex ,type)) *)
1343                  `(complex (- z (realpart w))
1344                            (- ,(coerce 0 ',type) (imagpart w))))
1345                ;; Multiply and divide two complex numbers.
1346                #!-complex-float-vops
1347                (deftransform * ((x y) ((complex ,type) (complex ,type)) *)
1348                  '(let* ((rx (realpart x))
1349                          (ix (imagpart x))
1350                          (ry (realpart y))
1351                          (iy (imagpart y)))
1352                     (complex (- (* rx ry) (* ix iy))
1353                              (+ (* rx iy) (* ix ry)))))
1354                (deftransform / ((x y) ((complex ,type) (complex ,type)) *)
1355                  #!-complex-float-vops
1356                  '(let* ((rx (realpart x))
1357                          (ix (imagpart x))
1358                          (ry (realpart y))
1359                          (iy (imagpart y)))
1360                     (if (> (abs ry) (abs iy))
1361                         (let* ((r (/ iy ry))
1362                                (dn (+ ry (* r iy))))
1363                           (complex (/ (+ rx (* ix r)) dn)
1364                                    (/ (- ix (* rx r)) dn)))
1365                         (let* ((r (/ ry iy))
1366                                (dn (+ iy (* r ry))))
1367                           (complex (/ (+ (* rx r) ix) dn)
1368                                    (/ (- (* ix r) rx) dn)))))
1369                  #!+complex-float-vops
1370                  `(let* ((cs (conjugate (sb!vm::swap-complex x)))
1371                          (ry (realpart y))
1372                          (iy (imagpart y)))
1373                     (if (> (abs ry) (abs iy))
1374                         (let* ((r (/ iy ry))
1375                                (dn (+ ry (* r iy))))
1376                           (/ (+ x (* cs r)) dn))
1377                         (let* ((r (/ ry iy))
1378                                (dn (+ iy (* r ry))))
1379                           (/ (+ (* x r) cs) dn)))))
1380                ;; Multiply a complex by a real or vice versa.
1381                #!-complex-float-vops
1382                (deftransform * ((w z) ((complex ,type) real) *)
1383                  '(complex (* (realpart w) z) (* (imagpart w) z)))
1384                #!-complex-float-vops
1385                (deftransform * ((z w) (real (complex ,type)) *)
1386                  '(complex (* (realpart w) z) (* (imagpart w) z)))
1387                ;; Divide a complex by a real or vice versa.
1388                #!-complex-float-vops
1389                (deftransform / ((w z) ((complex ,type) real) *)
1390                  '(complex (/ (realpart w) z) (/ (imagpart w) z)))
1391                (deftransform / ((x y) (,type (complex ,type)) *)
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 dn)
1399                                    (/ (- (* x r)) dn)))
1400                         (let* ((r (/ ry iy))
1401                                (dn (+ iy (* r ry))))
1402                           (complex (/ (* x r) dn)
1403                                    (/ (- x) dn)))))
1404                  #!+complex-float-vops
1405                  '(let* ((ry (realpart y))
1406                          (iy (imagpart y)))
1407                    (if (> (abs ry) (abs iy))
1408                        (let* ((r (/ iy ry))
1409                               (dn (+ ry (* r iy))))
1410                          (/ (complex x (- (* x r))) dn))
1411                        (let* ((r (/ ry iy))
1412                               (dn (+ iy (* r ry))))
1413                          (/ (complex (* x r) (- x)) dn)))))
1414                ;; conjugate of complex number
1415                #!-complex-float-vops
1416                (deftransform conjugate ((z) ((complex ,type)) *)
1417                  '(complex (realpart z) (- (imagpart z))))
1418                ;; CIS
1419                (deftransform cis ((z) ((,type)) *)
1420                  '(complex (cos z) (sin z)))
1421                ;; comparison
1422                #!-complex-float-vops
1423                (deftransform = ((w z) ((complex ,type) (complex ,type)) *)
1424                  '(and (= (realpart w) (realpart z))
1425                        (= (imagpart w) (imagpart z))))
1426                #!-complex-float-vops
1427                (deftransform = ((w z) ((complex ,type) real) *)
1428                  '(and (= (realpart w) z) (zerop (imagpart w))))
1429                #!-complex-float-vops
1430                (deftransform = ((w z) (real (complex ,type)) *)
1431                  '(and (= (realpart z) w) (zerop (imagpart z)))))))
1432
1433   (frob single-float)
1434   (frob double-float))
1435
1436 ;;; Here are simple optimizers for SIN, COS, and TAN. They do not
1437 ;;; produce a minimal range for the result; the result is the widest
1438 ;;; possible answer. This gets around the problem of doing range
1439 ;;; reduction correctly but still provides useful results when the
1440 ;;; inputs are union types.
1441 #-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
1442 (progn
1443 (defun trig-derive-type-aux (arg domain fun
1444                                  &optional def-lo def-hi (increasingp t))
1445   (etypecase arg
1446     (numeric-type
1447      (cond ((eq (numeric-type-complexp arg) :complex)
1448             (make-numeric-type :class (numeric-type-class arg)
1449                                :format (numeric-type-format arg)
1450                                :complexp :complex))
1451            ((numeric-type-real-p arg)
1452             (let* ((format (case (numeric-type-class arg)
1453                              ((integer rational) 'single-float)
1454                              (t (numeric-type-format arg))))
1455                    (bound-type (or format 'float)))
1456               ;; If the argument is a subset of the "principal" domain
1457               ;; of the function, we can compute the bounds because
1458               ;; the function is monotonic. We can't do this in
1459               ;; general for these periodic functions because we can't
1460               ;; (and don't want to) do the argument reduction in
1461               ;; exactly the same way as the functions themselves do
1462               ;; it.
1463               (if (csubtypep arg domain)
1464                   (let ((res-lo (bound-func fun (numeric-type-low arg) nil))
1465                         (res-hi (bound-func fun (numeric-type-high arg) nil)))
1466                     (unless increasingp
1467                       (rotatef res-lo res-hi))
1468                     (make-numeric-type
1469                      :class 'float
1470                      :format format
1471                      :low (coerce-numeric-bound res-lo bound-type)
1472                      :high (coerce-numeric-bound res-hi bound-type)))
1473                   (make-numeric-type
1474                    :class 'float
1475                    :format format
1476                    :low (and def-lo (coerce def-lo bound-type))
1477                    :high (and def-hi (coerce def-hi bound-type))))))
1478            (t
1479             (float-or-complex-float-type arg def-lo def-hi))))))
1480
1481 (defoptimizer (sin derive-type) ((num))
1482   (one-arg-derive-type
1483    num
1484    (lambda (arg)
1485      ;; Derive the bounds if the arg is in [-pi/2, pi/2].
1486      (trig-derive-type-aux
1487       arg
1488       (specifier-type `(float ,(- (/ pi 2)) ,(/ pi 2)))
1489       #'sin
1490       -1 1))
1491    #'sin))
1492
1493 (defoptimizer (cos derive-type) ((num))
1494   (one-arg-derive-type
1495    num
1496    (lambda (arg)
1497      ;; Derive the bounds if the arg is in [0, pi].
1498      (trig-derive-type-aux arg
1499                            (specifier-type `(float 0d0 ,pi))
1500                            #'cos
1501                            -1 1
1502                            nil))
1503    #'cos))
1504
1505 (defoptimizer (tan derive-type) ((num))
1506   (one-arg-derive-type
1507    num
1508    (lambda (arg)
1509      ;; Derive the bounds if the arg is in [-pi/2, pi/2].
1510      (trig-derive-type-aux arg
1511                            (specifier-type `(float ,(- (/ pi 2)) ,(/ pi 2)))
1512                            #'tan
1513                            nil nil))
1514    #'tan))
1515
1516 (defoptimizer (conjugate derive-type) ((num))
1517   (one-arg-derive-type num
1518     (lambda (arg)
1519       (flet ((most-negative-bound (l h)
1520                (and l h
1521                     (if (< (type-bound-number l) (- (type-bound-number h)))
1522                         l
1523                         (set-bound (- (type-bound-number h)) (consp h)))))
1524              (most-positive-bound (l h)
1525                (and l h
1526                     (if (> (type-bound-number h) (- (type-bound-number l)))
1527                         h
1528                         (set-bound (- (type-bound-number l)) (consp l))))))
1529         (if (numeric-type-real-p arg)
1530             (lvar-type num)
1531             (let ((low (numeric-type-low arg))
1532                   (high (numeric-type-high arg)))
1533               (let ((new-low (most-negative-bound low high))
1534                     (new-high (most-positive-bound low high)))
1535               (modified-numeric-type arg :low new-low :high new-high))))))
1536     #'conjugate))
1537
1538 (defoptimizer (cis derive-type) ((num))
1539   (one-arg-derive-type num
1540     (lambda (arg)
1541       (sb!c::specifier-type
1542        `(complex ,(or (numeric-type-format arg) 'float))))
1543     #'cis))
1544
1545 ) ; PROGN
1546 \f
1547 ;;;; TRUNCATE, FLOOR, CEILING, and ROUND
1548
1549 (macrolet ((define-frobs (fun ufun)
1550              `(progn
1551                 (defknown ,ufun (real) integer (movable foldable flushable))
1552                 (deftransform ,fun ((x &optional by)
1553                                     (* &optional
1554                                        (constant-arg (member 1))))
1555                   '(let ((res (,ufun x)))
1556                      (values res (- x res)))))))
1557   (define-frobs truncate %unary-truncate)
1558   (define-frobs round %unary-round))
1559
1560 (deftransform %unary-truncate ((x) (single-float))
1561   `(%unary-truncate/single-float x))
1562 (deftransform %unary-truncate ((x) (double-float))
1563   `(%unary-truncate/double-float x))
1564
1565 ;;; Convert (TRUNCATE x y) to the obvious implementation.
1566 ;;;
1567 ;;; ...plus hair: Insert explicit coercions to appropriate float types: Python
1568 ;;; is reluctant it generate explicit integer->float coercions due to
1569 ;;; precision issues (see SAFE-SINGLE-COERCION-P &co), but this is not an
1570 ;;; issue here as there is no DERIVE-TYPE optimizer on specialized versions of
1571 ;;; %UNARY-TRUNCATE, so the derived type of TRUNCATE remains the best we can
1572 ;;; do here -- which is fine. Also take care not to add unnecassary division
1573 ;;; or multiplication by 1, since we are not able to always eliminate them,
1574 ;;; depending on FLOAT-ACCURACY. Finally, leave out the secondary value when
1575 ;;; we know it is unused: COERCE is not flushable.
1576 (macrolet ((def (type other-float-arg-types)
1577              (let ((unary (symbolicate "%UNARY-TRUNCATE/" type))
1578                    (coerce (symbolicate "%" type)))
1579                `(deftransform truncate ((x &optional y)
1580                                         (,type
1581                                          &optional (or ,type ,@other-float-arg-types integer))
1582                                         * :result result)
1583                   (let* ((result-type (and result
1584                                            (lvar-derived-type result)))
1585                          (compute-all (and (values-type-p result-type)
1586                                            (not (type-single-value-p result-type)))))
1587                     (if (or (not y)
1588                             (and (constant-lvar-p y) (= 1 (lvar-value y))))
1589                         (if compute-all
1590                             `(let ((res (,',unary x)))
1591                                (values res (- x (,',coerce res))))
1592                             `(let ((res (,',unary x)))
1593                                ;; Dummy secondary value!
1594                                (values res x)))
1595                         (if compute-all
1596                             `(let* ((f (,',coerce y))
1597                                     (res (,',unary (/ x f))))
1598                                (values res (- x (* f (,',coerce res)))))
1599                             `(let* ((f (,',coerce y))
1600                                     (res (,',unary (/ x f))))
1601                                ;; Dummy secondary value!
1602                                (values res x)))))))))
1603   (def single-float ())
1604   (def double-float (single-float)))
1605
1606 (deftransform floor ((number &optional divisor)
1607                      (float &optional (or integer float)))
1608   (let ((defaulted-divisor (if divisor 'divisor 1)))
1609     `(multiple-value-bind (tru rem) (truncate number ,defaulted-divisor)
1610        (if (and (not (zerop rem))
1611                 (if (minusp ,defaulted-divisor)
1612                     (plusp number)
1613                     (minusp number)))
1614            (values (1- tru) (+ rem ,defaulted-divisor))
1615            (values tru rem)))))
1616
1617 (deftransform ceiling ((number &optional divisor)
1618                        (float &optional (or integer float)))
1619   (let ((defaulted-divisor (if divisor 'divisor 1)))
1620     `(multiple-value-bind (tru rem) (truncate number ,defaulted-divisor)
1621        (if (and (not (zerop rem))
1622                 (if (minusp ,defaulted-divisor)
1623                     (minusp number)
1624                     (plusp number)))
1625            (values (1+ tru) (- rem ,defaulted-divisor))
1626            (values tru rem)))))
1627
1628 (defknown %unary-ftruncate (real) float (movable foldable flushable))
1629 (defknown %unary-ftruncate/single (single-float) single-float
1630   (movable foldable flushable))
1631 (defknown %unary-ftruncate/double (double-float) double-float
1632   (movable foldable flushable))
1633
1634 (defun %unary-ftruncate/single (x)
1635   (declare (type single-float x))
1636   (declare (optimize speed (safety 0)))
1637   (let* ((bits (single-float-bits x))
1638          (exp (ldb sb!vm:single-float-exponent-byte bits))
1639          (biased (the single-float-exponent
1640                    (- exp sb!vm:single-float-bias))))
1641     (declare (type (signed-byte 32) bits))
1642     (cond
1643       ((= exp sb!vm:single-float-normal-exponent-max) x)
1644       ((<= biased 0) (* x 0f0))
1645       ((>= biased (float-digits x)) x)
1646       (t
1647        (let ((frac-bits (- (float-digits x) biased)))
1648          (setf bits (logandc2 bits (- (ash 1 frac-bits) 1)))
1649          (make-single-float bits))))))
1650
1651 (defun %unary-ftruncate/double (x)
1652   (declare (type double-float x))
1653   (declare (optimize speed (safety 0)))
1654   (let* ((high (double-float-high-bits x))
1655          (low (double-float-low-bits x))
1656          (exp (ldb sb!vm:double-float-exponent-byte high))
1657          (biased (the double-float-exponent
1658                    (- exp sb!vm:double-float-bias))))
1659     (declare (type (signed-byte 32) high)
1660              (type (unsigned-byte 32) low))
1661     (cond
1662       ((= exp sb!vm:double-float-normal-exponent-max) x)
1663       ((<= biased 0) (* x 0d0))
1664       ((>= biased (float-digits x)) x)
1665       (t
1666        (let ((frac-bits (- (float-digits x) biased)))
1667          (cond ((< frac-bits 32)
1668                 (setf low (logandc2 low (- (ash 1 frac-bits) 1))))
1669                (t
1670                 (setf low 0)
1671                 (setf high (logandc2 high (- (ash 1 (- frac-bits 32)) 1)))))
1672          (make-double-float high low))))))
1673
1674 (macrolet
1675     ((def (float-type fun)
1676          `(deftransform %unary-ftruncate ((x) (,float-type))
1677             '(,fun x))))
1678   (def single-float %unary-ftruncate/single)
1679   (def double-float %unary-ftruncate/double))