float-tran.lisp 72 KB
Newer Older
1
;;; -*- Mode: Lisp; Package: C; Log: code.log -*-
ram's avatar
ram committed
2 3
;;;
;;; **********************************************************************
ram's avatar
ram committed
4 5 6 7
;;; This code was written as part of the CMU Common Lisp project at
;;; Carnegie Mellon University, and has been placed in the public domain.
;;;
(ext:file-comment
8
  "$Header: src/compiler/float-tran.lisp $")
ram's avatar
ram committed
9
;;;
ram's avatar
ram committed
10 11 12 13 14 15 16 17
;;; **********************************************************************
;;;
;;; This file contains floating-point specific transforms, and may be somewhat
;;; implementation dependent in its assumptions of what the formats are.
;;;
;;; Author: Rob MacLachlan
;;; 
(in-package "C")
18
(intl:textdomain "cmucl")
ram's avatar
ram committed
19 20 21 22 23 24 25


;;;; Coercions:

(defknown %single-float (real) single-float (movable foldable flushable))
(defknown %double-float (real) double-float (movable foldable flushable))

26 27
(deftransform float ((n prototype) (* single-float) * :when :both)
  '(%single-float n))
ram's avatar
ram committed
28

29
(deftransform float ((n prototype) (* double-float) * :when :both)
ram's avatar
ram committed
30 31
  '(%double-float n))

32 33 34
(deftransform float ((n) *)
  `(if (floatp n) n (%single-float n)))

35
(deftransform %single-float ((n) (single-float) * :when :both)
ram's avatar
ram committed
36 37
  'n)

38
(deftransform %double-float ((n) (double-float) * :when :both)
ram's avatar
ram committed
39 40
  'n)

41 42 43 44
(defknown %complex-single-float (number) (complex single-float)
  (movable foldable flushable))
(defknown %complex-double-float (number) (complex double-float)
  (movable foldable flushable))
45
(defknown %complex-double-double-float (number) (complex double-double-float)
46 47 48 49 50 51 52 53
  (movable foldable flushable))

(macrolet
    ((frob (type)
       (let ((name (symbolicate "%COMPLEX-" type "-FLOAT"))
	     (convert (symbolicate "%" type "-FLOAT")))
	 `(progn
	    (defun ,name (n)
54 55
	      (declare (number n))
	      (etypecase n
56 57 58 59 60 61 62 63 64 65 66 67 68
		(real
		 (complex (,convert n)))
		(complex
		 (complex (,convert (realpart n))
			  (,convert (imagpart n))))))
	    (deftransform ,name ((n) ((complex ,(symbolicate type "-FLOAT"))) * :when :both)
	      'n)))))
  (frob single)
  (frob double)
  #+double-double
  (frob double-double))


ram's avatar
ram committed
69 70 71 72 73
(deftransform coerce ((n type) (* *) * :when :both)
  (unless (constant-continuation-p type)
    (give-up))
  `(the ,(continuation-value type)
	,(let ( (tspec (specifier-type (continuation-value type))) )
74 75 76 77
	   (cond #+double-double
		 ((csubtypep tspec (specifier-type 'double-double-float))
		  '(%double-double-float n))
		 ((csubtypep tspec (specifier-type 'double-float))	
ram's avatar
ram committed
78
		  '(%double-float n))	
79
		 ((csubtypep tspec (specifier-type 'single-float))
ram's avatar
ram committed
80
		  '(%single-float n))
81 82 83 84 85 86 87
		 #+double-double
		 ((csubtypep tspec (specifier-type '(complex double-double-float)))
		  '(%complex-double-double-float n))
		 ((csubtypep tspec (specifier-type '(complex double-float)))
		  '(%complex-double-float n))
		 ((csubtypep tspec (specifier-type '(complex single-float)))
		  '(%complex-single-float n))
ram's avatar
ram committed
88 89
		 (t
		  (give-up))))))
ram's avatar
ram committed
90

91
;;; Not strictly float functions, but primarily useful on floats:
92
;;;
93 94 95 96 97 98 99 100 101
(macrolet ((frob (fun ufun)
	     `(progn
		(defknown ,ufun (real) integer (movable foldable flushable))
		(deftransform ,fun ((x &optional by)
				    (* &optional
				       (constant-argument (member 1))))
		  '(let ((res (,ufun x)))
		     (values res (- x res)))))))
  (frob round %unary-round))
102

pw's avatar
pw committed
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
(defknown %unary-truncate (real) integer
	  (movable foldable flushable))

;; Convert (truncate x y) to the obvious implementation.  We only want
;; this when under certain conditions and let the generic truncate
;; handle the rest.  (Note: if y = 1, the divide and multiply by y
;; should be removed by other deftransforms.)

(deftransform truncate ((x &optional y)
			(float &optional (or float integer)))
  '(let ((res (%unary-truncate (/ x y))))
     (values res (- x (* y res)))))

(deftransform floor ((number &optional divisor)
		     (float &optional (or integer float)))
  '(multiple-value-bind (tru rem) (truncate number divisor)
    (if (and (not (zerop rem))
	     (if (minusp divisor)
		 (plusp number)
		 (minusp number)))
	(values (1- tru) (+ rem divisor))
	(values tru rem))))

(deftransform ceiling ((number &optional divisor)
		       (float &optional (or integer float)))
  '(multiple-value-bind (tru rem) (truncate number divisor)
    (if (and (not (zerop rem))
	     (if (minusp divisor)
		 (minusp number)
		 (plusp number)))
	(values (1+ tru) (- rem divisor))
	(values tru rem))))

136 137 138 139 140
(defknown %unary-ftruncate/single-float (single-float) single-float
	  (movable foldable flushable))
(defknown %unary-ftruncate/double-float (double-float) double-float
	  (movable foldable flushable))

141
(defknown %unary-ftruncate (real) float
142 143 144 145 146 147 148 149
	  (movable foldable flushable))

;; Convert (ftruncate x y) to the obvious implementation.  We only
;; want this under certain conditions and let the generic ftruncate
;; handle the rest.  (Note: if y = 1, the divide and multiply by y
;; should be removed by other deftransforms.)

(deftransform ftruncate ((x &optional (y 1))
150
			 (float &optional (or float integer)))
151 152 153
  '(let ((res (%unary-ftruncate (/ x y))))
     (values res (- x (* y res)))))

154
#+(or sparc (and x86 sse2))
155 156 157 158
(defknown fast-unary-ftruncate ((or single-float double-float))
  (or single-float double-float)
  (movable foldable flushable))

159
#+(or sparc (and x86 sse2))
160 161 162 163 164 165 166 167
(defoptimizer (fast-unary-ftruncate derive-type) ((f))
  (one-arg-derive-type f
		       #'(lambda (n)
			   (ftruncate-derive-type-quot-aux n
							   (specifier-type '(integer 1 1))
							   nil))
		       #'ftruncate))

168
;; Convert %unary-ftruncate to unary-ftruncate/{single,double}-float
169 170 171
;; if x is known to be of the right type.  Also, if the result is
;; known to fit in the same range as a (signed-byte 32), convert this
;; to %unary-truncate, which might be a single instruction, and float
172 173 174
;; the result.  However, for sparc and x86, we have a vop to do this
;; so call that, and for Sparc V9, we can actually handle a 64-bit
;; integer range.
175 176
(macrolet ((frob (ftype func)
	     `(deftransform %unary-ftruncate ((x) (,ftype))
177 178 179 180 181 182 183 184
	       (let* ((x-type (continuation-type x))
		      (lo (bound-value (numeric-type-low x-type)))
		      (hi (bound-value (numeric-type-high x-type)))
		      (limit-lo (- (ash 1 #-sparc-v9 31 #+sparc-v9 63)))
		      (limit-hi (ash 1 #-sparc-v9 31 #+sparc-v9 63)))
		 (if (and (numberp lo) (numberp hi)
			  (< limit-lo lo)
			  (< hi limit-hi))
185 186
		     #-(or sparc (and x86 sse2))
		     '(let ((result (coerce (%unary-truncate x) ',ftype)))
187 188 189
		       ;; Multiply by x when result is 0 so that we
		       ;; get the correct signed zero to match what
		       ;; ftruncate in float.lisp would return.
190 191 192 193 194 195 196 197
		       (if (zerop result)
			   (* result x)
			   result))
		     #+(or sparc (and x86 sse2))
		     '(let ((result (fast-unary-ftruncate x)))
		       (if (zerop result)
			   (* result x)
			   result))
198
		     '(,func x))))))
199 200
  (frob single-float %unary-ftruncate/single-float)
  (frob double-float %unary-ftruncate/double-float))
pw's avatar
pw committed
201

202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
;;; FROUND
#-x87
(progn
(deftransform fround ((x &optional (y 1))
		      ((or single-float double-float)
		       &optional (or single-float double-float integer)))
  '(let ((res (%unary-fround (/ x y))))
    (values res (- x (* y res)))))

(defknown %unary-fround (real) float
  (movable foldable flushable))

(defknown %unary-fround/single-float (single-float) single-float
  (movable foldable flushable))

(defknown %unary-fround/double-float (double-float) double-float
  (movable foldable flushable))

(deftransform %unary-fround ((x) (single-float))
  '(%unary-fround/single-float x))

(deftransform %unary-fround ((x) (double-float))
  '(%unary-fround/double-float x))

); not x87

ram's avatar
ram committed
228 229 230 231
;;; Random:
;;;
(macrolet ((frob (fun type)
	     `(deftransform random ((num &optional state)
232 233
				    (,type &optional *) *
				    :when :both)
ram's avatar
ram committed
234
		"use inline float operations"
235
		'(,fun num (or state *random-state*)))))
ram's avatar
ram committed
236 237 238
  (frob %random-single-float single-float)
  (frob %random-double-float double-float))

239
#-(or new-random random-mt19937 rand-xoroshiro)
ram's avatar
ram committed
240 241
(deftransform random ((num &optional state)
		      ((integer 1 #.random-fixnum-max) &optional *))
242
  _N"use inline fixnum operations"
ram's avatar
ram committed
243 244
  '(rem (random-chunk (or state *random-state*)) num))

245 246 247 248 249
;;; With the latest propagate-float-type code the compiler can inline
;;; truncate (signed-byte 32) allowing 31 bits, and (unsigned-byte 32)
;;; 32 bits on the x86. When not using the propagate-float-type
;;; feature the best size that can be inlined is 29 bits.  The choice
;;; shouldn't cause bootstrap problems just slow code.
250 251
#+new-random
(deftransform random ((num &optional state)
252
		      ((integer 1
pw's avatar
pw committed
253 254 255
				#+x86 #xffffffff
				#-x86 #x7fffffff
				)
256
		       &optional *))
257 258
  #+x86 (intl:gettext "use inline (unsigned-byte 32) operations")
  #-x86 (intl:gettext "use inline (signed-byte 32) operations")
259 260
  '(values (truncate (%random-double-float (coerce num 'double-float)
		      (or state *random-state*)))))
261

262
#+(or random-mt19937)
263 264
(deftransform random ((num &optional state)
		      ((integer 1 #.(expt 2 32)) &optional *))
265
  _N"use inline (unsigned-byte 32) operations"
266 267 268 269 270 271 272 273 274 275 276 277 278
  (let* ((num-type (continuation-type num))
	 (num-high (cond ((numeric-type-p num-type)
			  (numeric-type-high num-type))
			 ((union-type-p num-type)
			  ;; Find the maximum of the union type.  We
			  ;; know this works because if we're in this
			  ;; routine, NUM must be a subtype of
			  ;; (INTEGER 1 2^32), so each member of the
			  ;; union must be a subtype too.
			  (reduce #'max (union-type-types num-type)
				  :key #'numeric-type-high))
			 (t
			  (give-up)))))
rtoy's avatar
rtoy committed
279 280 281 282 283 284 285
    ;; Rather than doing (rem (random-chunk) num-high), we do,
    ;; essentially, (rem (* num-high (random-chunk)) #x100000000).  I
    ;; (rtoy) believe this approach doesn't have the bias issue with
    ;; doing rem.  This method works by treating (random-chunk) as if
    ;; it were a 32-bit fraction between 0 and 1, exclusive.  Multiply
    ;; this by num-high to get a random number between 0 and num-high,
    ;; This should have no bias.
286
    (cond ((constant-continuation-p num)
rtoy's avatar
rtoy committed
287 288 289 290 291
	   (if (= num-high (expt 2 32))
	       '(random-chunk (or state *random-state*))
	       '(values (bignum::%multiply 
			 (random-chunk (or state *random-state*))
			 num))))
292 293 294
	  ((< num-high (expt 2 32))
	   '(values (bignum::%multiply (random-chunk (or state *random-state*))
		     num)))
rtoy's avatar
rtoy committed
295 296 297 298 299
	  ((= num-high (expt 2 32))
	   '(if (= num (expt 2 32))
		(random-chunk (or state *random-state*))
		(values (bignum::%multiply (random-chunk (or state *random-state*))
					   num))))
300
	  (t
301
	   (error (intl:gettext "Shouldn't happen"))))))
302

ram's avatar
ram committed
303 304 305 306 307 308

;;;; Float accessors:

(defknown make-single-float ((signed-byte 32)) single-float
  (movable foldable flushable))

309
(defknown make-double-float ((signed-byte 32) (unsigned-byte 32)) double-float
ram's avatar
ram committed
310 311 312 313 314 315 316 317 318 319 320
  (movable foldable flushable))

(defknown single-float-bits (single-float) (signed-byte 32)
  (movable foldable flushable))

(defknown double-float-high-bits (double-float) (signed-byte 32)
  (movable foldable flushable))

(defknown double-float-low-bits (double-float) (unsigned-byte 32)
  (movable foldable flushable))

321
#+(or sparc ppc (and x86 sse2))
322 323 324 325
(defknown double-float-bits (double-float)
  (values (signed-byte 32) (unsigned-byte 32))
  (movable foldable flushable))

326 327 328 329 330 331 332 333 334 335 336 337 338 339 340
(deftransform float-sign ((float &optional float2)
			  (single-float &optional single-float) *)
  (if float2
      (let ((temp (gensym)))
	`(let ((,temp (abs float2)))
	  (if (minusp (single-float-bits float)) (- ,temp) ,temp)))
      '(if (minusp (single-float-bits float)) -1f0 1f0)))

(deftransform float-sign ((float &optional float2)
			  (double-float &optional double-float) *)
  (if float2
      (let ((temp (gensym)))
	`(let ((,temp (abs float2)))
	  (if (minusp (double-float-high-bits float)) (- ,temp) ,temp)))
      '(if (minusp (double-float-high-bits float)) -1d0 1d0)))
ram's avatar
ram committed
341 342 343 344


;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, SCALE-FLOAT:
;;;
345
;;;    Convert these operations to format specific versions when the format is
ram's avatar
ram committed
346 347 348
;;; known.
;;;

349
(deftype single-float-exponent ()
350 351 352
  `(integer ,(- vm:single-float-normal-exponent-min vm:single-float-bias
		vm:single-float-digits)
	    ,(- vm:single-float-normal-exponent-max vm:single-float-bias)))
353 354

(deftype double-float-exponent ()
355 356 357
  `(integer ,(- vm:double-float-normal-exponent-min vm:double-float-bias
		vm:double-float-digits)
	    ,(- vm:double-float-normal-exponent-max vm:double-float-bias)))
358 359 360


(deftype single-float-int-exponent ()
361 362 363 364
  `(integer ,(- vm:single-float-normal-exponent-min vm:single-float-bias
		(* vm:single-float-digits 2))
	    ,(- vm:single-float-normal-exponent-max vm:single-float-bias
		vm:single-float-digits)))
ram's avatar
ram committed
365 366

(deftype double-float-int-exponent ()
367 368 369 370
  `(integer ,(- vm:double-float-normal-exponent-min vm:double-float-bias
		(* vm:double-float-digits 2))
	    ,(- vm:double-float-normal-exponent-max vm:double-float-bias
		vm:double-float-digits)))
ram's avatar
ram committed
371 372

(deftype single-float-significand ()
373
  `(integer 0 (,(ash 1 vm:single-float-digits))))
ram's avatar
ram committed
374 375

(deftype double-float-significand ()
376
  `(integer 0 (,(ash 1 vm:double-float-digits))))
ram's avatar
ram committed
377 378

(defknown decode-single-float (single-float)
379 380 381
  (values (single-float 0.5f0 (1f0))
	  single-float-exponent
	  (member -1f0 1f0))
ram's avatar
ram committed
382 383 384
  (movable foldable flushable))

(defknown decode-double-float (double-float)
385 386 387
  (values (double-float 0.5d0 (1d0))
	  double-float-exponent
	  (member -1d0 1d0))
ram's avatar
ram committed
388 389 390 391 392 393 394 395
  (movable foldable flushable))

(defknown integer-decode-single-float (single-float)
  (values single-float-significand single-float-int-exponent (integer -1 1))
  (movable foldable flushable))

(defknown integer-decode-double-float (double-float)
  (values double-float-significand double-float-int-exponent (integer -1 1))
ram's avatar
ram committed
396
  (movable foldable flushable))
ram's avatar
ram committed
397 398 399 400 401 402 403

(defknown scale-single-float (single-float fixnum) single-float
  (movable foldable flushable))

(defknown scale-double-float (double-float fixnum) double-float
  (movable foldable flushable))

404
(deftransform decode-float ((x) (single-float) * :when :both)
ram's avatar
ram committed
405 406
  '(decode-single-float x))

407
(deftransform decode-float ((x) (double-float) * :when :both)
ram's avatar
ram committed
408 409
  '(decode-double-float x))

410
(deftransform integer-decode-float ((x) (single-float) * :when :both)
ram's avatar
ram committed
411 412
  '(integer-decode-single-float x))

413
(deftransform integer-decode-float ((x) (double-float) * :when :both)
ram's avatar
ram committed
414 415
  '(integer-decode-double-float x))

416
(deftransform scale-float ((f ex) (single-float *) * :when :both)
417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437
  (cond ((and (backend-featurep :x86)
	      (not (backend-featurep :sse2))
	      (csubtypep (continuation-type ex)
			 (specifier-type '(signed-byte 32)))
	      (not (byte-compiling)))
	 '(coerce (%scalbn (coerce f 'double-float) ex) 'single-float))
	((csubtypep (continuation-type ex)
		    (specifier-type `(integer #.(- vm:single-float-normal-exponent-min
						   vm:single-float-bias
						   vm:single-float-digits)
					      #.(- vm:single-float-normal-exponent-max
						   vm:single-float-bias
						   1))))
	 ;; The exponent is such that 2^ex will fit in a single-float.
	 ;; Thus, scale-float can be done multiplying by a suitable
	 ;; constant.
	 `(* f (kernel:make-single-float (dpb (+ ex (1+ vm:single-float-bias))
					      vm:single-float-exponent-byte
					      (kernel:single-float-bits 1f0)))))
	(t
	 '(scale-single-float f ex))))
ram's avatar
ram committed
438

439
(deftransform scale-float ((f ex) (double-float *) * :when :both)
440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460
  (cond ((and (backend-featurep :x86)
	      (not (backend-featurep :sse2))
	      (csubtypep (continuation-type ex)
			 (specifier-type '(signed-byte 32))))
	 '(%scalbn f ex))
	((csubtypep (continuation-type ex)
		    (specifier-type `(integer #.(- vm:double-float-normal-exponent-min
						   vm:double-float-bias
						   vm:double-float-digits)
					      #.(- vm:double-float-normal-exponent-max
						   vm:double-float-bias
						   1))))
	 ;; The exponent is such that 2^ex will fit in a double-float.
	 ;; Thus, scale-float can be done multiplying by a suitable
	 ;; constant.
	 `(* f (kernel:make-double-float (dpb (+ ex (1+ vm:double-float-bias))
					      vm:double-float-exponent-byte
					      (kernel::double-float-bits 1d0))
					 0)))
	(t
	 '(scale-double-float f ex))))
ram's avatar
ram committed
461 462 463 464 465 466

;;; toy@rtp.ericsson.se:
;;;
;;; Optimizers for scale-float.  If the float has bounds, new bounds
;;; are computed for the result, if possible.

467 468 469 470 471 472 473
(defun scale-float-derive-type-aux (f ex same-arg)
  (declare (ignore same-arg))
  (flet ((scale-bound (x n)
	   ;; We need to be a bit careful here and catch any overflows
	   ;; that might occur.  We can ignore underflows which become
	   ;; zeros.
	   (set-bound
474 475 476 477 478 479 480 481 482 483 484
	    (let ((value (handler-case
			     (scale-float (bound-value x) n)
			   (floating-point-overflow ()
			     nil))))
	      ;; This check is necessary for ppc because the current
	      ;; implementation on ppc doesn't signal floating-point
	      ;; overflow.  (How many other places do we need to check
	      ;; for this?)
	      (if (and (floatp value) (float-infinity-p value))
		  nil
		  value))
485 486 487 488 489 490 491 492 493 494 495 496
	    (consp x))))
    (when (and (numeric-type-p f) (numeric-type-p ex))
      (let ((f-lo (numeric-type-low f))
	    (f-hi (numeric-type-high f))
	    (ex-lo (numeric-type-low ex))
	    (ex-hi (numeric-type-high ex))
	    (new-lo nil)
	    (new-hi nil))
	(when (and f-hi ex-hi)
	  (setf new-hi (scale-bound f-hi ex-hi)))
	(when (and f-lo ex-lo)
	  (setf new-lo (scale-bound f-lo ex-lo)))
497 498 499 500 501 502 503
	;; We're computing bounds for scale-float.  Assume the bounds
	;; on f are fl and fh, and the bounds on ex are nl and nh.
	;; The resulting bound should be fl*2^nl and fh*2^nh.
	;; However, if fh is negative, and we get an underflow, we
	;; might get bounds like 0 and fh*2^nh < 0.  Our bounds are
	;; backwards.  Thus, swap the bounds to get the correct
	;; bounds.
504 505
	(when (and new-lo new-hi (< (bound-value new-hi)
				    (bound-value new-lo)))
506
	  (rotatef new-lo new-hi))
507 508 509 510 511 512
	(make-numeric-type :class (numeric-type-class f)
			   :format (numeric-type-format f)
			   :complexp :real
			   :low new-lo
			   :high new-hi)))))
;;;
513
(defoptimizer (scale-float derive-type) ((f ex))
514
  (two-arg-derive-type f ex #'scale-float-derive-type-aux
515
		       #'scale-float))
ram's avatar
ram committed
516
	     
pw's avatar
pw committed
517 518
;;; toy@rtp.ericsson.se:
;;;
ram's avatar
ram committed
519 520
;;; Defoptimizers for %single-float and %double-float.  This makes the
;;; FLOAT function return the correct ranges if the input has some
521
;;; defined range.  Quite useful if we want to convert some type of
ram's avatar
ram committed
522 523
;;; bounded integer into a float.

524 525
(macrolet
    ((frob (fun type)
526
       (let ((aux-name (symbolicate fun "-DERIVE-TYPE-AUX")))
527
	 `(progn
528
	   (defun ,aux-name (num)
529
	     ;; When converting a number to a float, the limits are
530
	     ;; the "same."  
531
	     (let* ((lo (bound-func #'(lambda (x)
532 533 534 535 536 537 538
					;; If we can't coerce it, we
					;; return a NIL for the bound.
					;; (Is IGNORE-ERRORS too
					;; heavy-handed?  Should we
					;; try to do something more
					;; fine-grained?)
					(ignore-errors (coerce x ',type)))
539 540
				    (numeric-type-low num)))
		    (hi (bound-func #'(lambda (x)
541
					(ignore-errors (coerce x ',type)))
542 543 544 545
				    (numeric-type-high num))))
	       (specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
	   
	   (defoptimizer (,fun derive-type) ((num))
dtc's avatar
dtc committed
546
	     (one-arg-derive-type num #',aux-name #',fun))))))
ram's avatar
ram committed
547 548
  (frob %single-float single-float)
  (frob %double-float double-float))
549

ram's avatar
ram committed
550

551 552 553 554
;;;; Float contagion:

;;; FLOAT-CONTAGION-ARG1, ARG2  --  Internal
;;;
555
;;;    Do some stuff to recognize when the loser is doing mixed float and
556 557
;;; rational arithmetic, or different float types, and fix it up.  If we don't,
;;; he won't even get so much as an efficency note.
558 559 560 561 562 563 564 565 566 567 568 569
;;;
(deftransform float-contagion-arg1 ((x y) * * :defun-only t :node node)
  `(,(continuation-function-name (basic-combination-fun node))
    (float x y) y))
;;;
(deftransform float-contagion-arg2 ((x y) * * :defun-only t :node node)
  `(,(continuation-function-name (basic-combination-fun node))
    x (float y x)))

(dolist (x '(+ * / -))
  (%deftransform x '(function (rational float) *) #'float-contagion-arg1)
  (%deftransform x '(function (float rational) *) #'float-contagion-arg2))
ram's avatar
ram committed
570

571 572 573 574 575 576 577 578 579 580 581
(dolist (x '(= < > + * / -))
  (%deftransform x '(function (single-float double-float) *)
		 #'float-contagion-arg1)
  (%deftransform x '(function (double-float single-float) *)
		 #'float-contagion-arg2))


;;; Prevent zerop, plusp, minusp from losing horribly.  We can't in general
;;; float rational args to comparison, since Common Lisp semantics says we are
;;; supposed to compare as rationals, but we can do it for any rational that
;;; has a precise representation as a float (such as 0).
ram's avatar
ram committed
582 583
;;;
(macrolet ((frob (op)
ram's avatar
ram committed
584
	     `(deftransform ,op ((x y) (float rational) * :when :both)
585
		(unless (constant-continuation-p y)
586
		  (give-up (intl:gettext "Can't open-code float to rational comparison.")))
587 588
		(let ((val (continuation-value y)))
		  (unless (eql (rational (float val)) val)
589
		    (give-up (intl:gettext "~S doesn't have a precise float representation.")
590 591
			     val)))
		`(,',op x (float y x)))))
ram's avatar
ram committed
592 593 594 595
  (frob <)
  (frob >)
  (frob =))

596 597 598 599 600 601
;; Convert (/ x n) to (* x (/ n)) when x is a float and n is a power
;; of two, because (/ n) can be reprsented exactly.
(deftransform / ((x y) (float float) * :when :both)
  (unless (constant-continuation-p y)
    (give-up))
  (let ((val (continuation-value y)))
Raymond Toy's avatar
Raymond Toy committed
602 603 604
    (unless (= (decode-float val) 0.5)
      (give-up))
    `(* x (float (/ ,val) x))))
605 606 607 608 609 610 611 612 613 614 615

;; Convert 2*x to x+x.
(deftransform * ((x y) (float real) * :when :both)
  (unless (constant-continuation-p y)
    (give-up))
  (let ((val (continuation-value y)))
    (unless (= val 2)
      (give-up))
    '(+ x x)))

	      
616 617

;;;; Irrational transforms:
618

rtoy's avatar
rtoy committed
619
(defknown (%tan %sinh %asinh %atanh %log %logb %log10 #+x87 %tan-quick)
620 621 622
	  (double-float) double-float
  (movable foldable flushable))

rtoy's avatar
rtoy committed
623
(defknown (%sin %cos %tanh #+x87 %sin-quick #+x87 %cos-quick)
ram's avatar
ram committed
624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641
    (double-float) (double-float -1.0d0 1.0d0)
    (movable foldable flushable))

(defknown (%asin %atan)
    (double-float) (double-float #.(- (/ pi 2)) #.(/ pi 2))
    (movable foldable flushable))
    
(defknown (%acos)
    (double-float) (double-float 0.0d0 #.pi)
    (movable foldable flushable))
    
(defknown (%cosh)
    (double-float) (double-float 1.0d0)
    (movable foldable flushable))

(defknown (%acosh %exp %sqrt)
    (double-float) (double-float 0.0d0)
    (movable foldable flushable))
642 643 644 645

(defknown %expm1
    (double-float) (double-float -1d0)
    (movable foldable flushable))
ram's avatar
ram committed
646 647 648

(defknown (%hypot)
    (double-float double-float) (double-float 0d0)
649 650
  (movable foldable flushable))

ram's avatar
ram committed
651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666
(defknown (%pow)
    (double-float double-float) double-float
  (movable foldable flushable))

(defknown (%atan2)
    (double-float double-float) (double-float #.(- pi) #.pi)
  (movable foldable flushable))

(defknown (%scalb)
    (double-float double-float) double-float
  (movable foldable flushable))

(defknown (%scalbn)
    (double-float (signed-byte 32)) double-float
    (movable foldable flushable))

667
(defknown (%log1p)
ram's avatar
ram committed
668 669 670
    (double-float) double-float
    (movable foldable flushable))

671 672 673
(dolist (stuff '((exp %exp *)
		 (log %log float)
		 (sqrt %sqrt float)
674
		 (sin %sin float)
675 676
		 (cos %cos float)
		 (tan %tan float)
677 678 679 680 681 682 683 684 685 686 687 688
		 (asin %asin float)
		 (acos %acos float)
		 (atan %atan *)
		 (sinh %sinh *)
		 (cosh %cosh *)
		 (tanh %tanh *)
		 (asinh %asinh *)
		 (acosh %acosh float)
		 (atanh %atanh float)))
  (destructuring-bind (name prim rtype) stuff
    (deftransform name ((x) '(single-float) rtype :eval-name t)
      `(coerce (,prim (coerce x 'double-float)) 'single-float))
689
    (deftransform name ((x) '(double-float) rtype :eval-name t :when :both)
690 691
      `(,prim x))))

Raymond Toy's avatar
Raymond Toy committed
692
(defknown (%sincos)
693 694 695 696 697
    (double-float) (values double-float double-float)
    (movable foldable flushable))

(deftransform cis ((x) (single-float) * :when :both)
  `(multiple-value-bind (s c)
Raymond Toy's avatar
Raymond Toy committed
698
       (%sincos (coerce x 'double-float))
699 700
     (complex (coerce c 'single-float)
	      (coerce s 'single-float))))
701 702

(deftransform cis ((x) (double-float) * :when :both)
703
  `(multiple-value-bind (s c)
Raymond Toy's avatar
Raymond Toy committed
704
       (%sincos x)
705
     (complex c s)))
706

ram's avatar
ram committed
707 708 709
;;; The argument range is limited on the x86 FP trig. functions. A
;;; post-test can detect a failure (and load a suitable result), but
;;; this test is avoided if possible.
710 711 712 713
;;;
;;; Simple tests show that sin/cos produce numbers greater than 1 when
;;; the arg >= 2^63.  tan produces floating-point invalid exceptions
;;; for arg >= 2^62.  So limit these to that range.
rtoy's avatar
rtoy committed
714
#+x87
715 716 717 718 719
(dolist (stuff '((sin %sin %sin-quick 63)
		 (cos %cos %cos-quick 63)
		 (tan %tan %tan-quick 62)))
  (destructuring-bind (name prim prim-quick limit)
      stuff
720
    (deftransform name ((x) '(single-float) '* :eval-name t)
721 722
      (if (and (backend-featurep :x86)
	       (not (backend-featurep :sse2)))
723
	  (cond ((csubtypep (continuation-type x)
724 725 726
			    (specifier-type `(single-float
					      (,(- (expt 2f0 limit)))
					      (,(expt 2f0 limit)))))
727 728 729 730
		 `(coerce (,prim-quick (coerce x 'double-float))
		   'single-float))
		(t 
		 (compiler-note
731
		  _N"Unable to avoid inline argument range check~@
732 733 734
                      because the argument range (~s) was not within 2^~D"
		  (type-specifier (continuation-type x))
		  limit)
735 736
		 `(coerce (,prim (coerce x 'double-float)) 'single-float)))
	  `(coerce (,prim (coerce x 'double-float)) 'single-float)))
ram's avatar
ram committed
737
    (deftransform name ((x) '(double-float) '* :eval-name t :when :both)
738 739
      (if (and (backend-featurep :x86)
	       (not (backend-featurep :sse2)))
740
	  (cond ((csubtypep (continuation-type x)
741 742 743
			    (specifier-type `(double-float
					      (,(- (expt 2d0 limit)))
					      (,(expt 2d0 limit)))))
744 745 746
		 `(,prim-quick x))
		(t 
		 (compiler-note
747
		  _N"Unable to avoid inline argument range check~@
748 749 750
                   because the argument range (~s) was not within 2^~D"
		  (type-specifier (continuation-type x))
		  limit)
751 752
		 `(,prim x)))
	  `(,prim x)))))
ram's avatar
ram committed
753

754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769
(deftransform atan ((x y) (single-float single-float) *)
  `(coerce (%atan2 (coerce x 'double-float) (coerce y 'double-float))
    'single-float))
(deftransform atan ((x y) (double-float double-float) * :when :both)
  `(%atan2 x y))

(deftransform expt ((x y) ((single-float 0f0) single-float) *)
  `(coerce (%pow (coerce x 'double-float) (coerce y 'double-float))
    'single-float))
(deftransform expt ((x y) ((double-float 0d0) double-float) * :when :both)
  `(%pow x y))
(deftransform expt ((x y) ((single-float 0f0) (signed-byte 32)) *)
  `(coerce (%pow (coerce x 'double-float) (coerce y 'double-float))
    'single-float))
(deftransform expt ((x y) ((double-float 0d0) (signed-byte 32)) * :when :both)
  `(%pow x (coerce y 'double-float)))
770

pw's avatar
pw committed
771 772 773 774
;;; ANSI says log with base zero returns zero.
(deftransform log ((x y) (float float) float)
  '(if (zerop y) y (/ (log x) (log y))))

775 776 777 778 779 780 781 782 783 784 785 786
(deftransform log ((x y) ((or (member 0f0) (single-float (0f0)))
			  (constant-argument number))
		   single-float)
  ;; Transform (log x 2) and (log x 10) to something simpler.
  (let ((y-val (continuation-value y)))
    (unless (and (not-more-contagious y x)
		 (or (= y-val 2)
		     (= y-val 10)))
      (give-up))
    (cond ((= y-val 10)
	   `(coerce (kernel:%log10 (float x 1d0)) 'single-float))
	  ((= y-val 2)
Raymond Toy's avatar
Raymond Toy committed
787
	   `(coerce (kernel:%log2 (float x 1d0)) 'single-float)))))
788 789 790 791 792 793 794 795 796 797 798 799 800

(deftransform log ((x y) ((or (member 0d0) (double-float 0d0))
			  (constant-argument number))
		   double-float)
  ;; Transform (log x 2) and (log x 10) to something simpler.
  (let ((y-val (continuation-value y)))
    (unless (and (not-more-contagious y x)
		 (or (= y-val 2)
		     (= y-val 10)))
      (give-up))
    (cond ((= y-val 10)
	   `(kernel:%log10 (float x 1d0)))
	  ((= y-val 2)
Raymond Toy's avatar
Raymond Toy committed
801
	   `(kernel:%log2 (float x 1d0))))))
ram's avatar
ram committed
802 803 804 805 806 807

;;; Handle some simple transformations
  
(deftransform abs ((x) ((complex double-float)) double-float :when :both)
  '(%hypot (realpart x) (imagpart x)))

808
(deftransform abs ((x) ((complex single-float)) single-float)
dtc's avatar
dtc committed
809 810
  '(coerce (%hypot (coerce (realpart x) 'double-float)
		   (coerce (imagpart x) 'double-float))
ram's avatar
ram committed
811 812
	  'single-float))

813 814 815 816 817 818 819 820 821
(deftransform abs ((x) (real) real)
  (let ((x-type (continuation-type x)))
    ;; If the arg is known to non-negative, we can just return the
    ;; arg.  However, (abs -0.0) is 0.0, so this transform only works
    ;; on floats that are known not to include negative zero.
    (if (csubtypep x-type (specifier-type '(or (rational 0) (float (0d0)) (member 0f0 0d0))))
	'x
	(give-up))))

ram's avatar
ram committed
822 823 824
(deftransform phase ((x) ((complex double-float)) double-float :when :both)
  '(%atan2 (imagpart x) (realpart x)))

825
(deftransform phase ((x) ((complex single-float)) single-float)
dtc's avatar
dtc committed
826 827
  '(coerce (%atan2 (coerce (imagpart x) 'double-float)
		   (coerce (realpart x) 'double-float))
ram's avatar
ram committed
828 829 830 831
	  'single-float))

(deftransform phase ((x) ((float)) float :when :both)
  '(if (minusp (float-sign x))
dtc's avatar
dtc committed
832 833
       (float pi x)
       (float 0 x)))
834 835

;;; The number is of type REAL.
836
(declaim (inline numeric-type-real-p))
837 838 839 840
(defun numeric-type-real-p (type)
  (and (numeric-type-p type)
       (eq (numeric-type-complexp type) :real)))

841 842 843
;;; Coerce a numeric type bound to the given type while handling
;;; exclusive bounds.
(defun coerce-numeric-bound (bound type)
844 845
  (when bound
    (if (consp bound)
846
	(list (coerce (car bound) type))
847
	(coerce bound type))))
848

849 850 851 852 853 854 855

;;;; Optimizers for elementary functions
;;;;
;;;; These optimizers compute the output range of the elementary
;;;; function, based on the domain of the input.
;;;;

856 857 858 859 860 861 862 863 864 865
;;; Generate a specifier for a complex type specialized to the same
;;; type as the argument.
(defun complex-float-type (arg)
  (declare (type numeric-type arg))
  (let* ((format (case (numeric-type-class arg)
		   ((integer rational) 'single-float)
		   (t (numeric-type-format arg))))
	 (float-type (or format 'float)))
    (specifier-type `(complex ,float-type))))

dtc's avatar
dtc committed
866 867 868
;;; Compute a specifier like '(or float (complex float)), except float
;;; should be the right kind of float.  Allow bounds for the float
;;; part too.
869
(defun float-or-complex-float-type (arg &optional lo hi)
870 871 872 873 874
  (declare (type numeric-type arg))
  (let* ((format (case (numeric-type-class arg)
		   ((integer rational) 'single-float)
		   (t (numeric-type-format arg))))
	 (float-type (or format 'float))
875 876
	 (lo (coerce-numeric-bound lo float-type))
	 (hi (coerce-numeric-bound hi float-type)))
877 878
    (specifier-type `(or (,float-type ,(or lo '*) ,(or hi '*))
		         (complex ,float-type)))))
dtc's avatar
dtc committed
879

880 881 882 883
;;; Domain-Subtype
;;;
;;; Test if the numeric-type ARG is within in domain specified by
;;; DOMAIN-LOW and DOMAIN-HIGH, consider negative and positive zero to
dtc's avatar
dtc committed
884 885 886 887 888
;;; be distinct as for the :negative-zero-is-not-zero feature. Note
;;; that only inclusive and open domain limits are handled as these
;;; are the only types of limits currently used. With the
;;; :negative-zero-is-not-zero feature this could be handled by the
;;; numeric subtype code in type.lisp.
889 890 891 892 893 894 895 896 897 898 899
;;;
(defun domain-subtypep (arg domain-low domain-high)
  (declare (type numeric-type arg)
	   (type (or real null) domain-low domain-high))
  (let* ((arg-lo (numeric-type-low arg))
	 (arg-lo-val (bound-value arg-lo))
	 (arg-hi (numeric-type-high arg))
	 (arg-hi-val (bound-value arg-hi)))
    ;; Check that the ARG bounds are correctly canonicalised.
    (when (and arg-lo (floatp arg-lo-val) (zerop arg-lo-val) (consp arg-lo)
	       (minusp (float-sign arg-lo-val)))
900
      (compiler-note _N"Float zero bound ~s not correctly canonicalised?" arg-lo)
dtc's avatar
dtc committed
901
      (setq arg-lo 0l0 arg-lo-val 0l0))
902 903
    (when (and arg-hi (zerop arg-hi-val) (floatp arg-hi-val) (consp arg-hi)
	       (plusp (float-sign arg-hi-val)))
904
      (compiler-note _N"Float zero bound ~s not correctly canonicalised?" arg-hi)
dtc's avatar
dtc committed
905 906 907 908 909 910 911 912 913 914 915 916 917
      (setq arg-hi -0l0 arg-hi-val -0l0))
    (flet ((fp-neg-zero-p (f)	; Is F -0.0?
	     (and (floatp f) (zerop f) (minusp (float-sign f))))
	   (fp-pos-zero-p (f)	; Is F +0.0? 
	     (and (floatp f) (zerop f) (plusp (float-sign f)))))
      (and (or (null domain-low)
	       (and arg-lo (>= arg-lo-val domain-low)
		    (not (and (fp-pos-zero-p domain-low)
			      (fp-neg-zero-p arg-lo)))))
	   (or (null domain-high)
	       (and arg-hi (<= arg-hi-val domain-high)
		    (not (and (fp-neg-zero-p domain-high)
			      (fp-pos-zero-p arg-hi)))))))))
ram's avatar
ram committed
918

919
;;; Elfun-Derive-Type-Simple
920
;;; 
921 922
;;; Handle monotonic functions of a single variable whose domain is
;;; possibly part of the real line.  ARG is the variable, FCN is the
923 924 925 926 927 928
;;; function, and DOMAIN is a specifier that gives the (real) domain
;;; of the function.  If ARG is a subset of the DOMAIN, we compute the
;;; bounds directly.  Otherwise, we compute the bounds for the
;;; intersection between ARG and DOMAIN, and then append a complex
;;; result, which occurs for the parts of ARG not in the DOMAIN.
;;;
929 930 931 932 933
;;; Negative and positive zero are considered distinct within
;;; DOMAIN-LOW and DOMAIN-HIGH, as for the :negative-zero-is-not-zero
;;; feature.
;;;
;;; DEFAULT-LOW and DEFAULT-HIGH are the lower and upper bounds if we
934 935
;;; can't compute the bounds using FCN.
;;;
936 937
(defun elfun-derive-type-simple (arg fcn domain-low domain-high
				     default-low default-high
938
				     &optional (increasingp t))
939
  (declare (type (or null real) domain-low domain-high))
940 941 942
  (etypecase arg
    (numeric-type
     (cond ((eq (numeric-type-complexp arg) :complex)
943
	    (complex-float-type arg))
944
	   ((numeric-type-real-p arg)
945 946 947 948 949 950 951 952
	    ;; The argument is real, so let's find the intersection
	    ;; between the argument and the domain of the function.
	    ;; We compute the bounds on the intersection, and for
	    ;; everything else, we return a complex number of the
	    ;; appropriate type.
	    (multiple-value-bind (intersection difference)
		(interval-intersection/difference
		 (numeric-type->interval arg)
953
		 (make-interval :low domain-low :high domain-high))
954 955 956 957 958 959
	      (cond
		(intersection
		 ;; Process the intersection.
		 (let* ((low (interval-low intersection))
			(high (interval-high intersection))
			(res-lo (or (bound-func fcn (if increasingp low high))
960
				    default-low))
961
			(res-hi (or (bound-func fcn (if increasingp high low))
962
				    default-high))
963 964 965 966 967 968 969 970 971 972 973 974 975 976
			;; Result specifier type.
			(format (case (numeric-type-class arg)
				  ((integer rational) 'single-float)
				  (t (numeric-type-format arg))))
			(bound-type (or format 'float))
			(result-type 
			 (make-numeric-type
			  :class 'float
			  :format format
			  :low (coerce-numeric-bound res-lo bound-type)
			  :high (coerce-numeric-bound res-hi bound-type))))
		   ;; If the ARG is a subset of the domain, we don't
		   ;; have to worry about the difference, because that
		   ;; can't occur.
977 978 979
		   (if (or (null difference)
			   ;; Check if the arg is within the domain.
			   (domain-subtypep arg domain-low domain-high))
980 981 982 983
		       result-type
		       (list result-type
			     (specifier-type `(complex ,bound-type))))))
		(t
984 985
		 ;; No intersection so the result must be purely complex.
		 (complex-float-type arg)))))
986
	   (t
987
	    (float-or-complex-float-type arg default-low default-high))))))
988 989

(macrolet
990 991
    ((frob (name domain-low domain-high def-low-bnd def-high-bnd
		 &key (increasingp t))
992
       (let ((num (gensym)))
993 994 995 996 997
	 `(defoptimizer (,name derive-type) ((,num))
	   (one-arg-derive-type
	    ,num
	    #'(lambda (arg)
		(elfun-derive-type-simple arg #',name
998 999
					  ,domain-low ,domain-high
					  ,def-low-bnd ,def-high-bnd
1000
					  ,increasingp))
dtc's avatar
dtc committed
1001
	    #',name)))))
ram's avatar
ram committed
1002 1003
  ;; These functions are easy because they are defined for the whole
  ;; real line.
1004 1005 1006 1007
  (frob exp nil nil 0 nil)
  (frob sinh nil nil nil nil)
  (frob tanh nil nil -1 1)
  (frob asinh nil nil nil nil)
ram's avatar
ram committed
1008 1009

  ;; These functions are only defined for part of the real line.  The
1010
  ;; condition selects the desired part of the line.  
1011
  (frob asin -1d0 1d0 (- (/ pi 2)) (/ pi 2))
1012 1013
  ;; Acos is monotonic decreasing, so we need to swap the function
  ;; values at the lower and upper bounds of the input domain.
1014 1015 1016 1017 1018 1019
  (frob acos -1d0 1d0 0 pi :increasingp nil)
  (frob acosh 1d0 nil nil nil)
  (frob atanh -1d0 1d0 -1 1)
  ;; Kahan says that (sqrt -0.0) is -0.0, so use a specifier that
  ;; includes -0.0.
  (frob sqrt -0d0 nil 0 nil))
1020
 
1021 1022 1023
;;; Compute bounds for (expt x y).  This should be easy since (expt x
;;; y) = (exp (* y (log x))).  However, computations done this way
;;; have too much roundoff.  Thus we have to do it the hard way.
1024
;;;  
1025 1026 1027 1028 1029 1030 1031 1032
(defun safe-expt (x y)
  (handler-case
      (expt x y)
    (error ()
      nil)))

;;; Handle the case when x >= 1
(defun interval-expt-> (x y)
1033
  (case (c::interval-range-info y 0d0)
1034 1035
    ('+
     ;; Y is positive and log X >= 0.  The range of exp(y * log(x)) is
dtc's avatar
dtc committed
1036 1037
     ;; obviously non-negative.  We just have to be careful for
     ;; infinite bounds (given by nil).
1038 1039 1040 1041
     (let ((lo (safe-expt (c::bound-value (c::interval-low x))
			  (c::bound-value (c::interval-low y))))
	   (hi (safe-expt (c::bound-value (c::interval-high x))
			  (c::bound-value (c::interval-high y)))))
dtc's avatar
dtc committed
1042
       (list (c::make-interval :low (or lo 1) :high hi))))
1043 1044
    ('-
     ;; Y is negative and log x >= 0.  The range of exp(y * log(x)) is
dtc's avatar
dtc committed
1045 1046
     ;; obviously [0, 1].  However, underflow (nil) means 0 is the
     ;; result
1047 1048 1049 1050
     (let ((lo (safe-expt (c::bound-value (c::interval-high x))
			  (c::bound-value (c::interval-low y))))
	   (hi (safe-expt (c::bound-value (c::interval-low x))
			  (c::bound-value (c::interval-high y)))))
dtc's avatar
dtc committed
1051
       (list (c::make-interval :low (or lo 0) :high (or hi 1)))))
1052 1053 1054 1055 1056 1057 1058
    (t
     ;; Split the interval in half
     (destructuring-bind (y- y+)
	 (c::interval-split 0 y t)
       (list (interval-expt-> x y-)
	     (interval-expt-> x y+))))))

1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113
;;; Handle the case when x < 0, and when y is known to be an integer.
;;; In this case, we can do something useful because the x^y is still
;;; a real number if x and y are.
(defun interval-expt-<-0 (x y)
  #+(or)
  (progn
    (format t "x = ~A~%" x)
    (format t "range-info y (~A) = ~A~%" y (interval-range-info y)))
  (flet ((handle-positive-power-0 (x y)
	   ;; -1 <= X <= 0 and Y is positive.  We need to consider if
	   ;; Y contains an odd integer or not.  Find the smallest
	   ;; even and odd integer (if possible) contained in Y.
	   (let* ((y-lo (bound-value (interval-low y)))
		  (min-odd (if (oddp y-lo)
			       y-lo
			       (let ((y-odd (1+ y-lo)))
				 (if (interval-contains-p y-odd y)
				     y-odd
				     nil))))
		  (min-even (if (evenp y-lo)
				y-lo
				(let ((y-even (1+ y-lo)))
				  (if (interval-contains-p y-even y)
				      y-even
				      nil)))))
	     (cond ((and min-odd min-even)
		    ;; The Y interval contains both even and odd
		    ;; integers.  Then the lower bound is (least
		    ;; x)^(least positive odd), because this
		    ;; creates the most negative value.  The upper
		    ;; is (most x)^(least positive even), because
		    ;; this is the most positive number.
		    ;;
		    ;; (Recall that if |x|<1, |x|^y gets smaller as y
		    ;; increases.)
		    (let ((lo (safe-expt (bound-value (interval-low x))
					 min-odd))
			  (hi (safe-expt (bound-value (interval-high x))
					 min-even)))
		      (list (make-interval :low lo :high hi))))
		   (min-odd
		    ;; Y consists of just one odd integer.
		    (assert (oddp min-odd))
		    (let ((lo (safe-expt (bound-value (interval-low x))
					 min-odd))
			  (hi (safe-expt (bound-value (interval-high x))
					 min-odd)))
		      (list (make-interval :low lo :high hi))))
		   (min-even
		    ;; Y consists of just one even integer.
		    (assert (evenp min-even))
		    (let ((lo (safe-expt (bound-value (interval-high x))
					 min-even))
			  (hi (safe-expt (bound-value (interval-low x))
					 min-even)))
1114 1115 1116 1117 1118
		      (list (make-interval :low lo :high hi))))
		   (t
		    ;; No mininum even or odd integer, so Y has no
		    ;; lower bound
		    (list (make-interval :low nil :high nil))))))
1119 1120 1121 1122
	 (handle-positive-power-1 (x y)
	   ;; X <= -1, Y is a positive integer.  Find the largest even
	   ;; and odd integer contained in Y, if possible.
	   (let* ((y-hi (bound-value (interval-high y)))
1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138
		  (max-odd (if y-hi
			       (if (oddp y-hi)
				   y-hi
				   (let ((y-odd (1- y-hi)))
				     (if (interval-contains-p y-odd y)
					 y-odd
					 nil)))
			       nil))
		  (max-even (if y-hi
				(if (evenp y-hi)
				    y-hi
				    (let ((y-even (1- y-hi)))
				      (if (interval-contains-p y-even y)
					  y-even
					  nil)))
				nil)))
1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167
	     ;; At least one of max-odd and max-even must be non-NIL!
	     (cond ((and max-odd max-even)
		    ;; The Y interval contains both even and odd
		    ;; integers.  Then the lower bound is (least
		    ;; x)^(most positive odd), because this
		    ;; creates the most negative value.  The upper
		    ;; is (least x)^(most positive even), because
		    ;; this is the most positive number.
		    ;;
		    (let ((lo (safe-expt (bound-value (interval-low x))
					 max-odd))
			  (hi (safe-expt (bound-value (interval-low x))
					 max-even)))
		      (list (make-interval :low lo :high hi))))
		   (max-odd
		    ;; Y consists of just one odd integer.
		    (assert (oddp max-odd))
		    (let ((lo (safe-expt (bound-value (interval-low x))
					 max-odd))
			  (hi (safe-expt (bound-value (interval-high x))
					 max-odd)))
		      (list (make-interval :low lo :high hi))))
		   (max-even
		    ;; Y consists of just one even integer.
		    (assert (evenp max-even))
		    (let ((lo (safe-expt (bound-value (interval-high x))
					 max-even))
			  (hi (safe-expt (bound-value (interval-low x))
					 max-even)))
1168 1169 1170 1171 1172
		      (list (make-interval :low lo :high hi))))
		   (t
		    ;; No maximum even or odd integer, which means y
		    ;; is no upper bound.
		    (list (make-interval :low nil :high nil)))))))
1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239
    ;; We need to split into x < -1 and -1 <= x <= 0, first.
    (case (interval-range-info x -1)
      ('+
       ;; -1 <= x <= 0
       #+(or)
       (format t "x range +~%")
       (case (interval-range-info y 0)