float-tran.lisp 91 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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#+double-double
(progn
(defknown %double-double-float (real)
  double-double-float
  (movable foldable flushable))

(deftransform float ((n prototype) (* double-double-float) * :when :both)
  '(%double-double-float n))

(deftransform %double-float ((n) (double-double-float) * :when :both)
  '(double-double-hi n))

(deftransform %single-float ((n) (double-double-float) * :when :both)
  '(float (double-double-hi n) 1f0))

(deftransform %double-double-float ((n) (double-double-float) * :when :both)
  'n)

#+nil
(defun %double-double-float (n)
  (make-double-double-float (float n 1d0) 0d0))

63
64
65
;; Moved to code/float.lisp, because we need this relatively early in
;; the build process to handle float and real types.
#+nil
66
67
68
(defun %double-double-float (n)
  (typecase n
    (fixnum
69
     (%make-double-double-float (float n 1d0) 0d0))
70
    (single-float
71
     (%make-double-double-float (float n 1d0) 0d0))
72
    (double-float
73
     (%make-double-double-float (float n 1d0) 0d0))
74
75
76
77
78
79
80
81
    (double-double-float
     n)
    (bignum
     (bignum:bignum-to-float n 'double-double-float))
    (ratio
     (kernel::float-ratio n 'double-double-float))))
); progn

82
83
84
85
(defknown %complex-single-float (number) (complex single-float)
  (movable foldable flushable))
(defknown %complex-double-float (number) (complex double-float)
  (movable foldable flushable))
86
(defknown %complex-double-double-float (number) (complex double-double-float)
87
88
89
90
91
92
93
94
  (movable foldable flushable))

(macrolet
    ((frob (type)
       (let ((name (symbolicate "%COMPLEX-" type "-FLOAT"))
	     (convert (symbolicate "%" type "-FLOAT")))
	 `(progn
	    (defun ,name (n)
95
96
	      (declare (number n))
	      (etypecase n
97
98
99
100
101
102
103
104
105
106
107
108
109
		(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
110
111
112
113
114
(deftransform coerce ((n type) (* *) * :when :both)
  (unless (constant-continuation-p type)
    (give-up))
  `(the ,(continuation-value type)
	,(let ( (tspec (specifier-type (continuation-value type))) )
115
116
117
118
	   (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
119
120
121
		  '(%double-float n))	
		 ((csubtypep tspec (specifier-type 'float))
		  '(%single-float n))
122
123
124
125
126
127
128
		 #+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
129
130
		 (t
		  (give-up))))))
ram's avatar
ram committed
131

132
;;; Not strictly float functions, but primarily useful on floats:
133
;;;
134
135
136
137
138
139
140
141
142
(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))
143

pw's avatar
pw committed
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
(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))))

177
178
179
180
181
(defknown %unary-ftruncate/single-float (single-float) single-float
	  (movable foldable flushable))
(defknown %unary-ftruncate/double-float (double-float) double-float
	  (movable foldable flushable))

182
(defknown %unary-ftruncate (real) float
183
184
185
186
187
188
189
190
	  (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))
191
			 (float &optional (or float integer)))
192
193
194
  '(let ((res (%unary-ftruncate (/ x y))))
     (values res (- x (* y res)))))

195
#+(or sparc (and x86 sse2))
196
197
198
199
(defknown fast-unary-ftruncate ((or single-float double-float))
  (or single-float double-float)
  (movable foldable flushable))

200
#+(or sparc (and x86 sse2))
201
202
203
204
205
206
207
208
(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))

209
;; Convert %unary-ftruncate to unary-ftruncate/{single,double}-float
210
211
212
;; 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
213
214
215
;; 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.
216
217
(macrolet ((frob (ftype func)
	     `(deftransform %unary-ftruncate ((x) (,ftype))
218
219
220
221
222
223
224
225
	       (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))
226
227
		     #-(or sparc (and x86 sse2))
		     '(let ((result (coerce (%unary-truncate x) ',ftype)))
228
229
230
		       ;; Multiply by x when result is 0 so that we
		       ;; get the correct signed zero to match what
		       ;; ftruncate in float.lisp would return.
231
232
233
234
235
236
237
238
		       (if (zerop result)
			   (* result x)
			   result))
		     #+(or sparc (and x86 sse2))
		     '(let ((result (fast-unary-ftruncate x)))
		       (if (zerop result)
			   (* result x)
			   result))
239
		     '(,func x))))))
240
241
  (frob single-float %unary-ftruncate/single-float)
  (frob double-float %unary-ftruncate/double-float))
pw's avatar
pw committed
242

243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
;;; 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
269
270
271
272
;;; Random:
;;;
(macrolet ((frob (fun type)
	     `(deftransform random ((num &optional state)
273
274
				    (,type &optional *) *
				    :when :both)
ram's avatar
ram committed
275
		"use inline float operations"
276
		'(,fun num (or state *random-state*)))))
ram's avatar
ram committed
277
278
279
  (frob %random-single-float single-float)
  (frob %random-double-float double-float))

280
#-(or new-random random-mt19937)
ram's avatar
ram committed
281
282
(deftransform random ((num &optional state)
		      ((integer 1 #.random-fixnum-max) &optional *))
283
  _N"use inline fixnum operations"
ram's avatar
ram committed
284
285
  '(rem (random-chunk (or state *random-state*)) num))

286
287
288
289
290
;;; 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.
291
292
#+new-random
(deftransform random ((num &optional state)
293
		      ((integer 1
pw's avatar
pw committed
294
295
296
				#+x86 #xffffffff
				#-x86 #x7fffffff
				)
297
		       &optional *))
298
299
  #+x86 (intl:gettext "use inline (unsigned-byte 32) operations")
  #-x86 (intl:gettext "use inline (signed-byte 32) operations")
300
301
  '(values (truncate (%random-double-float (coerce num 'double-float)
		      (or state *random-state*)))))
302
303
304
305

#+random-mt19937
(deftransform random ((num &optional state)
		      ((integer 1 #.(expt 2 32)) &optional *))
306
  _N"use inline (unsigned-byte 32) operations"
307
308
309
310
311
312
313
314
315
316
317
318
319
  (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
320
321
322
323
324
325
326
    ;; 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.
327
    (cond ((constant-continuation-p num)
rtoy's avatar
rtoy committed
328
329
330
331
332
	   (if (= num-high (expt 2 32))
	       '(random-chunk (or state *random-state*))
	       '(values (bignum::%multiply 
			 (random-chunk (or state *random-state*))
			 num))))
333
334
335
	  ((< num-high (expt 2 32))
	   '(values (bignum::%multiply (random-chunk (or state *random-state*))
		     num)))
rtoy's avatar
rtoy committed
336
337
338
339
340
	  ((= 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))))
341
	  (t
342
	   (error (intl:gettext "Shouldn't happen"))))))
343

ram's avatar
ram committed
344
345
346
347
348
349

;;;; Float accessors:

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

350
(defknown make-double-float ((signed-byte 32) (unsigned-byte 32)) double-float
ram's avatar
ram committed
351
352
353
354
355
356
357
358
359
360
361
  (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))

362
#+(or sparc ppc (and x86 sse2))
363
364
365
366
(defknown double-float-bits (double-float)
  (values (signed-byte 32) (unsigned-byte 32))
  (movable foldable flushable))

367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
#+double-double
(progn
(defknown double-double-float-p (t)
  boolean
  (movable foldable flushable))

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


(defknown double-double-hi (double-double-float)
  double-float
  (movable foldable flushable))

(defknown double-double-lo (double-double-float)
  double-float
  (movable foldable flushable))

) ; progn
387

388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
(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
403

404
#+double-double
405
406
407
408
409
410
411
412
413
414
415
(deftransform float-sign ((float &optional float2)
			  (double-double-float &optional double-double-float) *)
  (if float2
      (let ((temp (gensym)))
	`(let ((,temp (abs float2)))
	   (if (minusp (float-sign (double-double-hi float)))
	       (- ,temp)
	       ,temp)))
      '(if (minusp (float-sign (double-double-hi float))) -1w0 1w0)))

  
ram's avatar
ram committed
416
417
418

;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, SCALE-FLOAT:
;;;
419
;;;    Convert these operations to format specific versions when the format is
ram's avatar
ram committed
420
421
422
;;; known.
;;;

423
(deftype single-float-exponent ()
424
425
426
  `(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)))
427
428

(deftype double-float-exponent ()
429
430
431
  `(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)))
432
433
434


(deftype single-float-int-exponent ()
435
436
437
438
  `(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
439
440

(deftype double-float-int-exponent ()
441
442
443
444
  `(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
445
446

(deftype single-float-significand ()
447
  `(integer 0 (,(ash 1 vm:single-float-digits))))
ram's avatar
ram committed
448
449

(deftype double-float-significand ()
450
  `(integer 0 (,(ash 1 vm:double-float-digits))))
ram's avatar
ram committed
451
452

(defknown decode-single-float (single-float)
453
454
455
  (values (single-float 0.5f0 (1f0))
	  single-float-exponent
	  (member -1f0 1f0))
ram's avatar
ram committed
456
457
458
  (movable foldable flushable))

(defknown decode-double-float (double-float)
459
460
461
  (values (double-float 0.5d0 (1d0))
	  double-float-exponent
	  (member -1d0 1d0))
ram's avatar
ram committed
462
463
464
465
466
467
468
469
  (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
470
  (movable foldable flushable))
ram's avatar
ram committed
471
472
473
474
475
476
477

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

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

478
(deftransform decode-float ((x) (single-float) * :when :both)
ram's avatar
ram committed
479
480
  '(decode-single-float x))

481
(deftransform decode-float ((x) (double-float) * :when :both)
ram's avatar
ram committed
482
483
  '(decode-double-float x))

484
(deftransform integer-decode-float ((x) (single-float) * :when :both)
ram's avatar
ram committed
485
486
  '(integer-decode-single-float x))

487
(deftransform integer-decode-float ((x) (double-float) * :when :both)
ram's avatar
ram committed
488
489
  '(integer-decode-double-float x))

490
(deftransform scale-float ((f ex) (single-float *) * :when :both)
491
  (if (and (backend-featurep :x86)
492
	   (not (backend-featurep :sse2))
493
494
495
496
497
	   (csubtypep (continuation-type ex)
		      (specifier-type '(signed-byte 32)))
	   (not (byte-compiling)))
      '(coerce (%scalbn (coerce f 'double-float) ex) 'single-float)
      '(scale-single-float f ex)))
ram's avatar
ram committed
498

499
(deftransform scale-float ((f ex) (double-float *) * :when :both)
500
  (if (and (backend-featurep :x86)
501
	   (not (backend-featurep :sse2))
502
503
504
505
	   (csubtypep (continuation-type ex)
		      (specifier-type '(signed-byte 32))))
      '(%scalbn f ex)
      '(scale-double-float f ex)))
ram's avatar
ram committed
506
507
508
509
510
511

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

512
513
514
515
516
517
518
(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
519
520
521
522
523
524
525
526
527
528
529
	    (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))
530
531
532
533
534
535
536
537
538
539
540
541
	    (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)))
542
543
544
545
546
547
548
	;; 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.
549
550
	(when (and new-lo new-hi (< (bound-value new-hi)
				    (bound-value new-lo)))
551
	  (rotatef new-lo new-hi))
552
553
554
555
556
557
	(make-numeric-type :class (numeric-type-class f)
			   :format (numeric-type-format f)
			   :complexp :real
			   :low new-lo
			   :high new-hi)))))
;;;
558
(defoptimizer (scale-float derive-type) ((f ex))
559
  (two-arg-derive-type f ex #'scale-float-derive-type-aux
560
		       #'scale-float))
ram's avatar
ram committed
561
	     
pw's avatar
pw committed
562
563
;;; toy@rtp.ericsson.se:
;;;
ram's avatar
ram committed
564
565
;;; Defoptimizers for %single-float and %double-float.  This makes the
;;; FLOAT function return the correct ranges if the input has some
566
;;; defined range.  Quite useful if we want to convert some type of
ram's avatar
ram committed
567
568
;;; bounded integer into a float.

569
570
(macrolet
    ((frob (fun type)
571
       (let ((aux-name (symbolicate fun "-DERIVE-TYPE-AUX")))
572
	 `(progn
573
	   (defun ,aux-name (num)
574
	     ;; When converting a number to a float, the limits are
575
	     ;; the "same."  
576
	     (let* ((lo (bound-func #'(lambda (x)
577
578
579
580
581
582
583
					;; 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)))
584
585
				    (numeric-type-low num)))
		    (hi (bound-func #'(lambda (x)
586
					(ignore-errors (coerce x ',type)))
587
588
589
590
				    (numeric-type-high num))))
	       (specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
	   
	   (defoptimizer (,fun derive-type) ((num))
dtc's avatar
dtc committed
591
	     (one-arg-derive-type num #',aux-name #',fun))))))
ram's avatar
ram committed
592
593
  (frob %single-float single-float)
  (frob %double-float double-float))
594

ram's avatar
ram committed
595

596
597
598
599
;;;; Float contagion:

;;; FLOAT-CONTAGION-ARG1, ARG2  --  Internal
;;;
600
;;;    Do some stuff to recognize when the loser is doing mixed float and
601
602
;;; 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.
603
604
605
606
607
608
609
610
611
612
613
614
;;;
(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
615

616
617
618
619
620
621
622
623
624
625
626
(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
627
628
;;;
(macrolet ((frob (op)
ram's avatar
ram committed
629
	     `(deftransform ,op ((x y) (float rational) * :when :both)
630
		(unless (constant-continuation-p y)
631
		  (give-up (intl:gettext "Can't open-code float to rational comparison.")))
632
633
		(let ((val (continuation-value y)))
		  (unless (eql (rational (float val)) val)
634
		    (give-up (intl:gettext "~S doesn't have a precise float representation.")
635
636
			     val)))
		`(,',op x (float y x)))))
ram's avatar
ram committed
637
638
639
640
  (frob <)
  (frob >)
  (frob =))

641
642
643
644
645
646
;; 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
647
648
649
    (unless (= (decode-float val) 0.5)
      (give-up))
    `(* x (float (/ ,val) x))))
650
651
652
653
654
655
656
657
658
659
660

;; 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)))

	      
661
662

;;;; Irrational transforms:
663

rtoy's avatar
rtoy committed
664
(defknown (%tan %sinh %asinh %atanh %log %logb %log10 #+x87 %tan-quick)
665
666
667
	  (double-float) double-float
  (movable foldable flushable))

rtoy's avatar
rtoy committed
668
(defknown (%sin %cos %tanh #+x87 %sin-quick #+x87 %cos-quick)
ram's avatar
ram committed
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
    (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))
687
688
689
690

(defknown %expm1
    (double-float) (double-float -1d0)
    (movable foldable flushable))
ram's avatar
ram committed
691
692
693

(defknown (%hypot)
    (double-float double-float) (double-float 0d0)
694
695
  (movable foldable flushable))

ram's avatar
ram committed
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
(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))

712
(defknown (%log1p)
ram's avatar
ram committed
713
714
715
    (double-float) double-float
    (movable foldable flushable))

716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
(dolist (stuff '((exp %exp *)
		 (log %log float)
		 (sqrt %sqrt float)
		 (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))
731
    (deftransform name ((x) '(double-float) rtype :eval-name t :when :both)
732
733
      `(,prim x))))

734
#+(or (and linux x86) sparc)
735
736
737
738
739
740
741
742
(progn
(defknown (kernel::%sincos)
    (double-float) (values double-float double-float)
    (movable foldable flushable))

(deftransform cis ((x) (single-float) * :when :both)
  `(multiple-value-bind (s c)
       (kernel::%sincos (coerce x 'double-float))
743
744
     (complex (coerce c 'single-float)
	      (coerce s 'single-float))))
745
746

(deftransform cis ((x) (double-float) * :when :both)
747
748
749
  `(multiple-value-bind (s c)
       (kernel::%sincos x)
     (complex c s)))
750
751
752
753
754
755
756

#+double-double
(deftransform cis ((z) (double-double-float) *)
  ;; Cis.
  '(complex (cos z) (sin z)))
)

ram's avatar
ram committed
757
758
759
;;; 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.
760
761
762
763
;;;
;;; 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
764
#+x87
765
766
767
768
769
(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
770
    (deftransform name ((x) '(single-float) '* :eval-name t)
771
772
      (if (and (backend-featurep :x86)
	       (not (backend-featurep :sse2)))
773
	  (cond ((csubtypep (continuation-type x)
774
775
776
			    (specifier-type `(single-float
					      (,(- (expt 2f0 limit)))
					      (,(expt 2f0 limit)))))
777
778
779
780
		 `(coerce (,prim-quick (coerce x 'double-float))
		   'single-float))
		(t 
		 (compiler-note
781
		  _N"Unable to avoid inline argument range check~@
782
783
784
                      because the argument range (~s) was not within 2^~D"
		  (type-specifier (continuation-type x))
		  limit)
785
786
		 `(coerce (,prim (coerce x 'double-float)) 'single-float)))
	  `(coerce (,prim (coerce x 'double-float)) 'single-float)))
ram's avatar
ram committed
787
    (deftransform name ((x) '(double-float) '* :eval-name t :when :both)
788
789
      (if (and (backend-featurep :x86)
	       (not (backend-featurep :sse2)))
790
	  (cond ((csubtypep (continuation-type x)
791
792
793
			    (specifier-type `(double-float
					      (,(- (expt 2d0 limit)))
					      (,(expt 2d0 limit)))))
794
795
796
		 `(,prim-quick x))
		(t 
		 (compiler-note
797
		  _N"Unable to avoid inline argument range check~@
798
799
800
                   because the argument range (~s) was not within 2^~D"
		  (type-specifier (continuation-type x))
		  limit)
801
802
		 `(,prim x)))
	  `(,prim x)))))
ram's avatar
ram committed
803

804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
(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)))
820

pw's avatar
pw committed
821
822
823
824
;;; ANSI says log with base zero returns zero.
(deftransform log ((x y) (float float) float)
  '(if (zerop y) y (/ (log x) (log y))))

ram's avatar
ram committed
825
826
827
828
829
830

;;; Handle some simple transformations
  
(deftransform abs ((x) ((complex double-float)) double-float :when :both)
  '(%hypot (realpart x) (imagpart x)))

831
(deftransform abs ((x) ((complex single-float)) single-float)
dtc's avatar
dtc committed
832
833
  '(coerce (%hypot (coerce (realpart x) 'double-float)
		   (coerce (imagpart x) 'double-float))
ram's avatar
ram committed
834
835
	  'single-float))

836
837
838
839
840
841
842
843
844
(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
845
846
847
(deftransform phase ((x) ((complex double-float)) double-float :when :both)
  '(%atan2 (imagpart x) (realpart x)))

848
(deftransform phase ((x) ((complex single-float)) single-float)
dtc's avatar
dtc committed
849
850
  '(coerce (%atan2 (coerce (imagpart x) 'double-float)
		   (coerce (realpart x) 'double-float))
ram's avatar
ram committed
851
852
853
854
	  'single-float))

(deftransform phase ((x) ((float)) float :when :both)
  '(if (minusp (float-sign x))
dtc's avatar
dtc committed
855
856
       (float pi x)
       (float 0 x)))
857
858

;;; The number is of type REAL.
859
(declaim (inline numeric-type-real-p))
860
861
862
863
(defun numeric-type-real-p (type)
  (and (numeric-type-p type)
       (eq (numeric-type-complexp type) :real)))

864
865
866
;;; Coerce a numeric type bound to the given type while handling
;;; exclusive bounds.
(defun coerce-numeric-bound (bound type)
867
868
  (when bound
    (if (consp bound)
869
	(list (coerce (car bound) type))
870
	(coerce bound type))))
871

872
873
874
875
876
877
878

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

879
880
881
882
883
884
885
886
887
888
;;; 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
889
890
891
;;; Compute a specifier like '(or float (complex float)), except float
;;; should be the right kind of float.  Allow bounds for the float
;;; part too.
892
(defun float-or-complex-float-type (arg &optional lo hi)
893
894
895
896
897
  (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))
898
899
	 (lo (coerce-numeric-bound lo float-type))
	 (hi (coerce-numeric-bound hi float-type)))
900
901
    (specifier-type `(or (,float-type ,(or lo '*) ,(or hi '*))
		         (complex ,float-type)))))
dtc's avatar
dtc committed
902

903
904
905
906
;;; 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
907
908
909
910
911
;;; 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.
912
913
914
915
916
917
918
919
920
921
922
;;;
(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)))
923
      (compiler-note _N"Float zero bound ~s not correctly canonicalised?" arg-lo)
dtc's avatar
dtc committed
924
      (setq arg-lo 0l0 arg-lo-val 0l0))
925
926
    (when (and arg-hi (zerop arg-hi-val) (floatp arg-hi-val) (consp arg-hi)
	       (plusp (float-sign arg-hi-val)))
927
      (compiler-note _N"Float zero bound ~s not correctly canonicalised?" arg-hi)
dtc's avatar
dtc committed
928
929
930
931
932
933
934
935
936
937
938
939
940
      (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
941

942
;;; Elfun-Derive-Type-Simple
943
;;; 
944
945
;;; Handle monotonic functions of a single variable whose domain is
;;; possibly part of the real line.  ARG is the variable, FCN is the
946
947
948
949
950
951
;;; 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.
;;;
952
953
954
955
956
;;; 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
957
958
;;; can't compute the bounds using FCN.
;;;
959
960
(defun elfun-derive-type-simple (arg fcn domain-low domain-high
				     default-low default-high
961
				     &optional (increasingp t))
962
  (declare (type (or null real) domain-low domain-high))
963
964
965
  (etypecase arg
    (numeric-type
     (cond ((eq (numeric-type-complexp arg) :complex)
966
	    (complex-float-type arg))
967
	   ((numeric-type-real-p arg)
968
969
970
971
972
973
974
975
	    ;; 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)
976
		 (make-interval :low domain-low :high domain-high))
977
978
979
980
981
982
	      (cond
		(intersection
		 ;; Process the intersection.
		 (let* ((low (interval-low intersection))
			(high (interval-high intersection))
			(res-lo (or (bound-func fcn (if increasingp low high))
983
				    default-low))
984
			(res-hi (or (bound-func fcn (if increasingp high low))
985
				    default-high))
986
987
988
989
990
991
992
993
994
995
996
997
998
999
			;; 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.
1000
1001
1002
		   (if (or (null difference)
			   ;; Check if the arg is within the domain.
			   (domain-subtypep arg domain-low domain-high))
1003
1004
1005
1006
		       result-type
		       (list result-type
			     (specifier-type `(complex ,bound-type))))))
		(t
1007
1008
		 ;; No intersection so the result must be purely complex.
		 (complex-float-type arg)))))
1009
	   (t
1010
	    (float-or-complex-float-type arg default-low default-high))))))
1011
1012

(macrolet
1013
1014
    ((frob (name domain-low domain-high def-low-bnd def-high-bnd
		 &key (increasingp t))
1015
       (let ((num (gensym)))
1016
1017
1018
1019
1020
	 `(defoptimizer (,name derive-type) ((,num))
	   (one-arg-derive-type
	    ,num
	    #'(lambda (arg)
		(elfun-derive-type-simple arg #',name
1021
1022
					  ,domain-low ,domain-high
					  ,def-low-bnd ,def-high-bnd
1023
					  ,increasingp))
dtc's avatar
dtc committed
1024
	    #',name)))))
ram's avatar
ram committed
1025
1026
  ;; These functions are easy because they are defined for the whole
  ;; real line.
1027
1028
1029
1030
  (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
1031
1032

  ;; These functions are only defined for part of the real line.  The
1033
  ;; condition selects the desired part of the line.  
1034
  (frob asin -1d0 1d0 (- (/ pi 2)) (/ pi 2))
1035
1036
  ;; Acos is monotonic decreasing, so we need to swap the function
  ;; values at the lower and upper bounds of the input domain.
1037
1038
1039
1040
1041
1042
  (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))
1043
 
1044
1045
1046
;;; 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.
1047
;;;  
1048
1049
1050
1051
1052
1053
1054
1055
(defun safe-expt (x y)
  (handler-case
      (expt x y)
    (error ()
      nil)))

;;; Handle the case when x >= 1
(defun interval-expt-> (x y)
1056
  (case (c::interval-range-info y 0d0)
1057
1058
    ('+
     ;; Y is positive and log X >= 0.  The range of exp(y * log(x)) is
dtc's avatar
dtc committed
1059
1060
     ;; obviously non-negative.  We just have to be careful for
     ;; infinite bounds (given by nil).
1061
1062
1063
1064
     (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
1065
       (list (c::make-interval :low (or lo 1) :high hi))))
1066
1067
    ('-
     ;; Y is negative and log x >= 0.  The range of exp(y * log(x)) is
dtc's avatar
dtc committed
1068
1069
     ;; obviously [0, 1].  However, underflow (nil) means 0 is the
     ;; result
1070
1071
1072
1073
     (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
1074
       (list (c::make-interval :low (or lo 0) :high (or hi 1)))))
1075
1076
1077
1078
1079
1080
1081
    (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+))))))

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
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
;;; 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)))
1137
1138
1139
1140
1141
		      (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))))))
1142
1143
1144
1145
	 (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)))
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
		  (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)))
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
	     ;; 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)))
1191
1192
1193
1194
1195
		      (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)))))))
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
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
    ;; 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)
	 ('+
	  (handle-positive-power-0 x y))
	 ('-
	  ;; Y is negative.  We should do something better
	  ;; than this because there's an extra rounding which
	  ;; we shouldn't do.
	  #+(or)
	  (format t "Handle y neg~%")
	  (let ((unit (make-interval :low 1 :high 1))
		(result (handle-positive-power-0 x (interval-neg y))))
	    #+(or)
	    (format t "result = ~A~%" result)
	    (mapcar #'(lambda (r)
			(interval-div unit r))
		    result)))
	 (t
	  ;; Split the interval and try again.  Since we know y is an
	  ;; integer, we don't need interval-split.  Also we want to
	  ;; handle an exponent of 0 ourselves as a special case.
	  (multiple-value-bind (y- y+)
	      (values (make-interval :low (interval-low y)
				     :high -1)
		      (make-interval :low 1
				     :high (interval-high y)))
	    (append (list (make-interval :low 1 :high 1))
		    (interval-expt-<-0 x y-)
		    (interval-expt-<-0 x y+))))))
      ('-
       ;; x < -1
       (case (c::interval-range-info y)
	 ('+
	  ;; Y is positive.  We need to consider if Y contains an
	  ;; odd integer or not.
	  ;;
	  (handle-positive-power-1 x y))
	 ('-
	  ;; Y is negative.  Do this in a better way
	  (let ((unit (make-interval :low 1 :high 1))
		(result (handle-positive-power-1 x (interval-neg y))))
	    (mapcar #'(lambda (r)
			(interval-div unit r))
		    result)))
	 (t
	  ;; Split the interval and try again.
	  #+(or)
	  (format t "split y ~A~%" y)
	  (multiple-value-bind (y- y+)
	      (values (make-interval :low (interval-low y) :high -1)
		      (make-interval :low 1 :high (interval-high y)))
	    (append (list (make-interval :low 1 :high 1))
		    (interval-expt-<-0 x y-)
		    (interval-expt-<-0 x y+))))))
      (t
       #+(or)
       (format t "splitting x ~A~%" x)
       (destructuring-bind (neg pos)
	   (interval-split -1 x t t)
	 (append (interval-expt-<-0 neg y)
		 (interval-expt-<-0 pos y)))))))

dtc's avatar
dtc committed
1263
;;; Handle the case when x <= 1
1264
(defun interval-expt-< (x y &optional integer-power-p)
1265
  (case (c::interval-range-info x 0d0)
1266
    ('+
dtc's avatar
dtc committed
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
     ;; The case of 0 <= x <= 1 is easy
     (case (c::interval-range-info y)
       ('+
	;; Y is positive and log X <= 0.  The range of exp(y * log(x)) is
	;; obviously [0, 1].  We just have to be careful for infinite bounds
	;; (given by nil).
	(let ((lo (safe-expt (c::bound-value (c::interval-low x))
			     (c::bound-value (c::interval-high y))))
	      (hi (safe-expt (c::bound-value (c::interval-high x))
			     (c::bound-value (c::interval-low y)))))
1277
1278
1279
1280
1281
1282
1283
1284
1285
	  ;; If the low bound, LO, is NIL, that means the we have
	  ;; +0.0^inf, which is +0.0, but NIL is returned by
	  ;; SAFE-EXPT.  That means the result is includes +0.0.  Make
	  ;; it so by returning a member type and an exclusive
	  ;; interval.
	  (if lo
	      (list (c::make-interval :low lo :high (or hi 1)))
	      (list (c::make-interval :low (list 0) :high (or hi 1))
		    (c::make-member-type :members (list 0))))))
dtc's avatar
dtc committed
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
       ('-
	;; Y is negative and log x <= 0.  The range of exp(y * log(x)) is
	;; obviously [1, inf].
	(let ((hi (safe-expt (c::bound-value (c::interval-low x))
			     (c::bound-value (c::interval-low y))))
	      (lo (safe-expt (c::bound-value (c::interval-high x))
			     (c::bound-value (c::interval-high y)))))
	  (list (c::make-interval :low (or lo 1) :high hi))))
       (t
	;; Split the interval in half
	(destructuring-bind (y- y+)
	    (c::interval-split 0 y t)
1298
1299
	  (list (interval-expt-< x y- integer-power-p)
		(interval-expt-< x y+ integer-power-p))))))
1300
    ('-
1301
1302
     ;; The case where x <= 0.
     (cond (integer-power-p
1303
1304
1305
1306
1307
	    (interval-expt-<-0 x y))
	   (t
	    ;; Y is not an integer.  Just give up and return an
	    ;; unbounded interval.
	    (list (c::make-interval :low nil :high nil)))))
1308
    (t
dtc's avatar
dtc committed
1309
1310
     (destructuring-bind (neg pos)
	 (interval-split 0 x t t)
1311
1312
       (append (interval-expt-< neg y integer-power-p)
	       (interval-expt-< pos y integer-power-p))))))
1313
1314

;;; Compute bounds for (expt x y)
dtc's avatar
dtc committed
1315

1316
(defun interval-expt (x y &optional integer-power-p)