float-tran.lisp 71.9 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: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/float-tran.lisp,v 1.114 2007/06/21 16:22:24 rtoy Exp $")
ram's avatar
ram committed
9
;;;
ram's avatar
ram committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
;;; **********************************************************************
;;;
;;; 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")


;;;; Coercions:

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

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

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

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

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

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

40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#+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))

(defun %double-double-float (n)
  (typecase n
    (fixnum
65
     (%make-double-double-float (float n 1d0) 0d0))
66
    (single-float
67
     (%make-double-double-float (float n 1d0) 0d0))
68
    (double-float
69
     (%make-double-double-float (float n 1d0) 0d0))
70
71
72
73
74
75
76
77
    (double-double-float
     n)
    (bignum
     (bignum:bignum-to-float n 'double-double-float))
    (ratio
     (kernel::float-ratio n 'double-double-float))))
); progn

ram's avatar
ram committed
78
79
80
81
82
(deftransform coerce ((n type) (* *) * :when :both)
  (unless (constant-continuation-p type)
    (give-up))
  `(the ,(continuation-value type)
	,(let ( (tspec (specifier-type (continuation-value type))) )
83
84
85
86
	   (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
87
88
89
90
91
		  '(%double-float n))	
		 ((csubtypep tspec (specifier-type 'float))
		  '(%single-float n))
		 (t
		  (give-up))))))
ram's avatar
ram committed
92

93
;;; Not strictly float functions, but primarily useful on floats:
94
;;;
95
96
97
98
99
100
101
102
103
(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))
104

pw's avatar
pw committed
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
136
137
(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))))

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

143
(defknown %unary-ftruncate (real) float
144
145
146
147
148
149
150
151
	  (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))
152
			 (float &optional (or float integer)))
153
154
155
  '(let ((res (%unary-ftruncate (/ x y))))
     (values res (- x (* y res)))))

156
157
158
159
160
#+sparc
(defknown fast-unary-ftruncate ((or single-float double-float))
  (or single-float double-float)
  (movable foldable flushable))

161
162
163
164
165
166
167
168
169
#+sparc
(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))

170
;; Convert %unary-ftruncate to unary-ftruncate/{single,double}-float
171
172
173
174
175
176
;; 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
;; the result.  However, for sparc, we have a vop to do this so call
;; that, and for Sparc V9, we can actually handle a 64-bit integer
;; range.
177
178
179

(macrolet ((frob (ftype func)
	     `(deftransform %unary-ftruncate ((x) (,ftype))
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
	       (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))
		     #-sparc '(let ((result (coerce (%unary-truncate x) ',ftype)))
			        (if (zerop result)
				    (* result x)
				    result))
		     #+sparc '(let ((result (fast-unary-ftruncate x)))
			        (if (zerop result)
				    (* result x)
				    result))
		     '(,func x))))))
197
198
  (frob single-float %unary-ftruncate/single-float)
  (frob double-float %unary-ftruncate/double-float))
pw's avatar
pw committed
199

ram's avatar
ram committed
200
201
202
203
;;; Random:
;;;
(macrolet ((frob (fun type)
	     `(deftransform random ((num &optional state)
204
205
				    (,type &optional *) *
				    :when :both)
ram's avatar
ram committed
206
		"use inline float operations"
207
		'(,fun num (or state *random-state*)))))
ram's avatar
ram committed
208
209
210
  (frob %random-single-float single-float)
  (frob %random-double-float double-float))

211
#-(or new-random random-mt19937)
ram's avatar
ram committed
212
213
214
215
216
(deftransform random ((num &optional state)
		      ((integer 1 #.random-fixnum-max) &optional *))
  "use inline fixnum operations"
  '(rem (random-chunk (or state *random-state*)) num))

217
218
219
220
221
;;; 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.
222
223
#+new-random
(deftransform random ((num &optional state)
224
		      ((integer 1
pw's avatar
pw committed
225
226
227
				#+x86 #xffffffff
				#-x86 #x7fffffff
				)
228
229
230
		       &optional *))
  #+x86 "use inline (unsigned-byte 32) operations"
  #-x86 "use inline (signed-byte 32) operations"
231
232
  '(values (truncate (%random-double-float (coerce num 'double-float)
		      (or state *random-state*)))))
233
234
235
236
237

#+random-mt19937
(deftransform random ((num &optional state)
		      ((integer 1 #.(expt 2 32)) &optional *))
  "use inline (unsigned-byte 32) operations"
238
239
240
241
242
243
244
245
246
247
248
249
250
  (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)))))
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
    (cond ((constant-continuation-p num)
	   ;; Check the worst case sum abs error for the random number
	   ;; expectations.
	   (let ((rem (rem (expt 2 32) num-high)))
	     (unless (< (/ (* 2 rem (- num-high rem)) num-high (expt 2 32))
			(expt 2 (- kernel::random-integer-extra-bits)))
	       (give-up "The random number expectations are inaccurate."))
	     (if (= num-high (expt 2 32))
		 '(random-chunk (or state *random-state*))
		 #-x86 '(rem (random-chunk (or state *random-state*)) num)
		 #+x86
		 ;; Use multiplication which is faster.
		 '(values (bignum::%multiply 
			   (random-chunk (or state *random-state*))
			   num)))))
	  ((> num-high random-fixnum-max)
	   (give-up "The range is too large to assure an accurate result."))
	  #+x86
	  ((< num-high (expt 2 32))
	   '(values (bignum::%multiply (random-chunk (or state *random-state*))
		     num)))
	  (t
	   '(rem (random-chunk (or state *random-state*)) num)))))

ram's avatar
ram committed
275
276
277
278
279
280

;;;; Float accessors:

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

281
(defknown make-double-float ((signed-byte 32) (unsigned-byte 32)) double-float
ram's avatar
ram committed
282
283
284
285
286
287
288
289
290
291
292
  (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))

293
#+(or sparc ppc)
294
295
296
297
(defknown double-float-bits (double-float)
  (values (signed-byte 32) (unsigned-byte 32))
  (movable foldable flushable))

298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
#+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
318

319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
(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
334

335
336
337
338
339
340
341
342
343
344
345
(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
346
347
348

;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, SCALE-FLOAT:
;;;
349
;;;    Convert these operations to format specific versions when the format is
ram's avatar
ram committed
350
351
352
;;; known.
;;;

353
(deftype single-float-exponent ()
354
355
356
  `(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)))
357
358

(deftype double-float-exponent ()
359
360
361
  `(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)))
362
363
364


(deftype single-float-int-exponent ()
365
366
367
368
  `(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
369
370

(deftype double-float-int-exponent ()
371
372
373
374
  `(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
375
376

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

(deftype double-float-significand ()
380
  `(integer 0 (,(ash 1 vm:double-float-digits))))
ram's avatar
ram committed
381
382

(defknown decode-single-float (single-float)
383
384
385
  (values (single-float 0.5f0 (1f0))
	  single-float-exponent
	  (member -1f0 1f0))
ram's avatar
ram committed
386
387
388
  (movable foldable flushable))

(defknown decode-double-float (double-float)
389
390
391
  (values (double-float 0.5d0 (1d0))
	  double-float-exponent
	  (member -1d0 1d0))
ram's avatar
ram committed
392
393
394
395
396
397
398
399
  (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
400
  (movable foldable flushable))
ram's avatar
ram committed
401
402
403
404
405
406
407

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

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

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

411
(deftransform decode-float ((x) (double-float) * :when :both)
ram's avatar
ram committed
412
413
  '(decode-double-float x))

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

417
(deftransform integer-decode-float ((x) (double-float) * :when :both)
ram's avatar
ram committed
418
419
  '(integer-decode-double-float x))

420
(deftransform scale-float ((f ex) (single-float *) * :when :both)
421
422
423
424
425
426
  (if (and (backend-featurep :x86)
	   (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
427

428
(deftransform scale-float ((f ex) (double-float *) * :when :both)
429
430
431
432
433
  (if (and (backend-featurep :x86)
	   (csubtypep (continuation-type ex)
		      (specifier-type '(signed-byte 32))))
      '(%scalbn f ex)
      '(scale-double-float f ex)))
ram's avatar
ram committed
434
435
436
437
438
439

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

440
441
442
443
444
445
446
(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
447
448
449
450
451
452
453
454
455
456
457
	    (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))
458
459
460
461
462
463
464
465
466
467
468
469
	    (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)))
470
471
472
473
474
475
476
	;; 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.
477
478
	(when (and new-lo new-hi (< (bound-value new-hi)
				    (bound-value new-lo)))
479
	  (rotatef new-lo new-hi))
480
481
482
483
484
485
	(make-numeric-type :class (numeric-type-class f)
			   :format (numeric-type-format f)
			   :complexp :real
			   :low new-lo
			   :high new-hi)))))
;;;
486
(defoptimizer (scale-float derive-type) ((f ex))
487
  (two-arg-derive-type f ex #'scale-float-derive-type-aux
488
		       #'scale-float t))
ram's avatar
ram committed
489
	     
pw's avatar
pw committed
490
491
;;; toy@rtp.ericsson.se:
;;;
ram's avatar
ram committed
492
493
;;; Defoptimizers for %single-float and %double-float.  This makes the
;;; FLOAT function return the correct ranges if the input has some
494
;;; defined range.  Quite useful if we want to convert some type of
ram's avatar
ram committed
495
496
;;; bounded integer into a float.

497
498
(macrolet
    ((frob (fun type)
499
       (let ((aux-name (symbolicate fun "-DERIVE-TYPE-AUX")))
500
	 `(progn
501
	   (defun ,aux-name (num)
502
	     ;; When converting a number to a float, the limits are
503
	     ;; the "same."  
504
	     (let* ((lo (bound-func #'(lambda (x)
505
506
507
508
509
510
511
					;; 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)))
512
513
				    (numeric-type-low num)))
		    (hi (bound-func #'(lambda (x)
514
					(ignore-errors (coerce x ',type)))
515
516
517
518
				    (numeric-type-high num))))
	       (specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
	   
	   (defoptimizer (,fun derive-type) ((num))
dtc's avatar
dtc committed
519
	     (one-arg-derive-type num #',aux-name #',fun))))))
ram's avatar
ram committed
520
521
  (frob %single-float single-float)
  (frob %double-float double-float))
522

ram's avatar
ram committed
523

524
525
526
527
;;;; Float contagion:

;;; FLOAT-CONTAGION-ARG1, ARG2  --  Internal
;;;
528
;;;    Do some stuff to recognize when the loser is doing mixed float and
529
530
;;; 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.
531
532
533
534
535
536
537
538
539
540
541
542
;;;
(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
543

544
545
546
547
548
549
550
551
552
553
554
(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
555
556
;;;
(macrolet ((frob (op)
ram's avatar
ram committed
557
	     `(deftransform ,op ((x y) (float rational) * :when :both)
558
559
		(unless (constant-continuation-p y)
		  (give-up "Can't open-code float to rational comparison."))
560
561
562
563
564
		(let ((val (continuation-value y)))
		  (unless (eql (rational (float val)) val)
		    (give-up "~S doesn't have a precise float representation."
			     val)))
		`(,',op x (float y x)))))
ram's avatar
ram committed
565
566
567
568
  (frob <)
  (frob >)
  (frob =))

569
570

;;;; Irrational transforms:
571

572
(defknown (%tan %sinh %asinh %atanh %log %logb %log10 %tan-quick)
573
574
575
	  (double-float) double-float
  (movable foldable flushable))

576
(defknown (%sin %cos %tanh %sin-quick %cos-quick)
ram's avatar
ram committed
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
    (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))
595
596
597
598

(defknown %expm1
    (double-float) (double-float -1d0)
    (movable foldable flushable))
ram's avatar
ram committed
599
600
601

(defknown (%hypot)
    (double-float double-float) (double-float 0d0)
602
603
  (movable foldable flushable))

ram's avatar
ram committed
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
(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))

620
(defknown (%log1p)
ram's avatar
ram committed
621
622
623
    (double-float) double-float
    (movable foldable flushable))

624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
(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))
639
    (deftransform name ((x) '(double-float) rtype :eval-name t :when :both)
640
641
      `(,prim x))))

ram's avatar
ram committed
642
643
644
645
646
647
648
;;; 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.
(dolist (stuff '((sin %sin %sin-quick)
		 (cos %cos %cos-quick)
		 (tan %tan %tan-quick)))
  (destructuring-bind (name prim prim-quick) stuff
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
    (deftransform name ((x) '(single-float) '* :eval-name t)
      (if (backend-featurep :x86)
	  (cond ((csubtypep (continuation-type x)
			    (specifier-type '(single-float
					      (#.(- (expt 2f0 64)))
					      (#.(expt 2f0 64)))))
		 `(coerce (,prim-quick (coerce x 'double-float))
		   'single-float))
		(t 
		 (compiler-note
		  "Unable to avoid inline argument range check~@
                      because the argument range (~s) was not within 2^64"
		  (type-specifier (continuation-type x)))
		 `(coerce (,prim (coerce x 'double-float)) 'single-float)))
	  `(coerce (,prim (coerce x 'double-float)) 'single-float)))
ram's avatar
ram committed
664
    (deftransform name ((x) '(double-float) '* :eval-name t :when :both)
665
666
667
668
669
670
671
672
673
674
675
676
677
      (if (backend-featurep :x86)
	  (cond ((csubtypep (continuation-type x)
			    (specifier-type '(double-float
					      (#.(- (expt 2d0 64)))
					      (#.(expt 2d0 64)))))
		 `(,prim-quick x))
		(t 
		 (compiler-note
		  "Unable to avoid inline argument range check~@
                   because the argument range (~s) was not within 2^64"
		  (type-specifier (continuation-type x)))
		 `(,prim x)))
	  `(,prim x)))))
ram's avatar
ram committed
678

679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
(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)))
695

pw's avatar
pw committed
696
697
698
699
;;; 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
700
701
702
703
704
705

;;; Handle some simple transformations
  
(deftransform abs ((x) ((complex double-float)) double-float :when :both)
  '(%hypot (realpart x) (imagpart x)))

706
(deftransform abs ((x) ((complex single-float)) single-float)
dtc's avatar
dtc committed
707
708
  '(coerce (%hypot (coerce (realpart x) 'double-float)
		   (coerce (imagpart x) 'double-float))
ram's avatar
ram committed
709
710
	  'single-float))

711
712
713
714
715
716
717
718
719
(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
720
721
722
(deftransform phase ((x) ((complex double-float)) double-float :when :both)
  '(%atan2 (imagpart x) (realpart x)))

723
(deftransform phase ((x) ((complex single-float)) single-float)
dtc's avatar
dtc committed
724
725
  '(coerce (%atan2 (coerce (imagpart x) 'double-float)
		   (coerce (realpart x) 'double-float))
ram's avatar
ram committed
726
727
728
729
	  'single-float))

(deftransform phase ((x) ((float)) float :when :both)
  '(if (minusp (float-sign x))
dtc's avatar
dtc committed
730
731
       (float pi x)
       (float 0 x)))
732
733

;;; The number is of type REAL.
734
(declaim (inline numeric-type-real-p))
735
736
737
738
(defun numeric-type-real-p (type)
  (and (numeric-type-p type)
       (eq (numeric-type-complexp type) :real)))

739
740
741
;;; Coerce a numeric type bound to the given type while handling
;;; exclusive bounds.
(defun coerce-numeric-bound (bound type)
742
743
  (when bound
    (if (consp bound)
744
	(list (coerce (car bound) type))
745
	(coerce bound type))))
746

747
748
749
750
751
752
753

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

754
755
756
757
758
759
760
761
762
763
;;; 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
764
765
766
;;; Compute a specifier like '(or float (complex float)), except float
;;; should be the right kind of float.  Allow bounds for the float
;;; part too.
767
(defun float-or-complex-float-type (arg &optional lo hi)
768
769
770
771
772
  (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))
773
774
	 (lo (coerce-numeric-bound lo float-type))
	 (hi (coerce-numeric-bound hi float-type)))
775
776
    (specifier-type `(or (,float-type ,(or lo '*) ,(or hi '*))
		         (complex ,float-type)))))
dtc's avatar
dtc committed
777

778
779
780
781
;;; 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
782
783
784
785
786
;;; 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.
787
788
789
790
791
792
793
794
795
796
797
798
;;;
(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)))
      (compiler-note "Float zero bound ~s not correctly canonicalised?" arg-lo)
dtc's avatar
dtc committed
799
      (setq arg-lo 0l0 arg-lo-val 0l0))
800
801
802
    (when (and arg-hi (zerop arg-hi-val) (floatp arg-hi-val) (consp arg-hi)
	       (plusp (float-sign arg-hi-val)))
      (compiler-note "Float zero bound ~s not correctly canonicalised?" arg-hi)
dtc's avatar
dtc committed
803
804
805
806
807
808
809
810
811
812
813
814
815
      (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
816

817
;;; Elfun-Derive-Type-Simple
818
;;; 
819
820
;;; Handle monotonic functions of a single variable whose domain is
;;; possibly part of the real line.  ARG is the variable, FCN is the
821
822
823
824
825
826
;;; 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.
;;;
827
828
829
830
831
;;; 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
832
833
;;; can't compute the bounds using FCN.
;;;
834
835
(defun elfun-derive-type-simple (arg fcn domain-low domain-high
				     default-low default-high
836
				     &optional (increasingp t))
837
  (declare (type (or null real) domain-low domain-high))
838
839
840
  (etypecase arg
    (numeric-type
     (cond ((eq (numeric-type-complexp arg) :complex)
841
	    (complex-float-type arg))
842
	   ((numeric-type-real-p arg)
843
844
845
846
847
848
849
850
	    ;; 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)
851
		 (make-interval :low domain-low :high domain-high))
852
853
854
855
856
857
	      (cond
		(intersection
		 ;; Process the intersection.
		 (let* ((low (interval-low intersection))
			(high (interval-high intersection))
			(res-lo (or (bound-func fcn (if increasingp low high))
858
				    default-low))
859
			(res-hi (or (bound-func fcn (if increasingp high low))
860
				    default-high))
861
862
863
864
865
866
867
868
869
870
871
872
873
874
			;; 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.
875
876
877
		   (if (or (null difference)
			   ;; Check if the arg is within the domain.
			   (domain-subtypep arg domain-low domain-high))
878
879
880
881
		       result-type
		       (list result-type
			     (specifier-type `(complex ,bound-type))))))
		(t
882
883
		 ;; No intersection so the result must be purely complex.
		 (complex-float-type arg)))))
884
	   (t
885
	    (float-or-complex-float-type arg default-low default-high))))))
886
887

(macrolet
888
889
    ((frob (name domain-low domain-high def-low-bnd def-high-bnd
		 &key (increasingp t))
890
       (let ((num (gensym)))
891
892
893
894
895
	 `(defoptimizer (,name derive-type) ((,num))
	   (one-arg-derive-type
	    ,num
	    #'(lambda (arg)
		(elfun-derive-type-simple arg #',name
896
897
					  ,domain-low ,domain-high
					  ,def-low-bnd ,def-high-bnd
898
					  ,increasingp))
dtc's avatar
dtc committed
899
	    #',name)))))
ram's avatar
ram committed
900
901
  ;; These functions are easy because they are defined for the whole
  ;; real line.
902
903
904
905
  (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
906
907

  ;; These functions are only defined for part of the real line.  The
908
  ;; condition selects the desired part of the line.  
909
  (frob asin -1d0 1d0 (- (/ pi 2)) (/ pi 2))
910
911
  ;; Acos is monotonic decreasing, so we need to swap the function
  ;; values at the lower and upper bounds of the input domain.
912
913
914
915
916
917
  (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))
918
 
919
920
921
;;; 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.
922
;;;  
923
924
925
926
927
928
929
930
(defun safe-expt (x y)
  (handler-case
      (expt x y)
    (error ()
      nil)))

;;; Handle the case when x >= 1
(defun interval-expt-> (x y)
931
  (case (c::interval-range-info y 0d0)
932
933
    ('+
     ;; Y is positive and log X >= 0.  The range of exp(y * log(x)) is
dtc's avatar
dtc committed
934
935
     ;; obviously non-negative.  We just have to be careful for
     ;; infinite bounds (given by nil).
936
937
938
939
     (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
940
       (list (c::make-interval :low (or lo 1) :high hi))))
941
942
    ('-
     ;; Y is negative and log x >= 0.  The range of exp(y * log(x)) is
dtc's avatar
dtc committed
943
944
     ;; obviously [0, 1].  However, underflow (nil) means 0 is the
     ;; result
945
946
947
948
     (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
949
       (list (c::make-interval :low (or lo 0) :high (or hi 1)))))
950
951
952
953
954
955
956
    (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+))))))

dtc's avatar
dtc committed
957
;;; Handle the case when x <= 1
958
(defun interval-expt-< (x y)
959
  (case (c::interval-range-info x 0d0)
960
    ('+
dtc's avatar
dtc committed
961
962
963
964
965
966
967
968
969
970
     ;; 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)))))
971
972
973
974
975
976
977
978
979
	  ;; 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
980
981
982
983
984
985
986
987
988
989
990
991
992
993
       ('-
	;; 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)
	  (list (interval-expt-< x y-)
		(interval-expt-< x y+))))))
994
    ('-
dtc's avatar
dtc committed
995
996
997
998
     ;; The case where x <= 0.  Y MUST be an INTEGER for this to
     ;; work!  The calling function must insure this! For now we'll
     ;; just return the appropriate unbounded float type.
     (list (c::make-interval :low nil :high nil)))
999
    (t
dtc's avatar
dtc committed
1000
1001
1002
1003
     (destructuring-bind (neg pos)
	 (interval-split 0 x t t)
       (list (interval-expt-< neg y)
	     (interval-expt-< pos y))))))
1004
1005

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

1007
(defun interval-expt (x y)
dtc's avatar
dtc committed
1008
1009
1010
  (case (interval-range-info x 1)
    ('+
     ;; X >= 1
1011
	 (interval-expt-> x y))
dtc's avatar
dtc committed
1012
1013
1014
1015
    ('-
     ;; X <= 1
     (interval-expt-< x y))
    (t
1016
     (destructuring-bind (left right)
dtc's avatar
dtc committed
1017
	 (interval-split 1 x t t)
1018
1019
       (list (interval-expt left y)
	     (interval-expt right y))))))
1020
1021

(defun fixup-interval-expt (bnd x-int y-int x-type y-type)
dtc's avatar
dtc committed
1022
  (declare (ignore x-int))
1023
1024
  ;; Figure out what the return type should be, given the argument
  ;; types and bounds and the result type and bounds.
rtoy's avatar
rtoy committed
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
  (flet ((low-bnd (b)
	   (etypecase b
	     (member-type
	      (reduce #'min (member-type-members b)))
	     (interval
	      (interval-low b))))
	 (hi-bnd (b)
	   (etypecase b
	     (member-type
	      (reduce #'max (member-type-members b)))
	     (interval
	      (interval-high b)))))
    (cond ((csubtypep x-type (specifier-type 'integer))
	   ;; An integer to some power.  Cases to consider:
	   (case (numeric-type-class y-type)
	     (integer
	      ;; Positive integer to an integer power is either an
	      ;; integer or a rational.
	      (let ((lo (or (low-bnd bnd) '*))
		    (hi (or (hi-bnd bnd) '*)))
		(if (and (interval-low y-int)
			 (>= (bound-value (interval-low y-int)) 0))
		    (specifier-type `(integer ,lo ,hi))
		    (specifier-type `(rational ,lo ,hi)))))
	     (rational
	      ;; Positive integer to rational power is either a rational
	      ;; or a single-float.
	      (let* ((lo (low-bnd bnd))
		     (hi (hi-bnd bnd))
		     (int-lo (if lo
				 (floor (bound-value lo))
				 '*))
		     (int-hi (if hi
				 (ceiling (bound-value hi))
				 '*))
		     (f-lo (if lo
			       (bound-func #'float lo)
1062
			       '*))
rtoy's avatar
rtoy committed
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
		     (f-hi (if hi
			       (bound-func #'float hi)
			       '*)))
		(specifier-type `(or (rational ,int-lo ,int-hi)
				     (single-float ,f-lo, f-hi)))))
	     (float
	      ;; Positive integer to a float power is a float
	      (let ((res (copy-numeric-type y-type)))
		(setf (numeric-type-low res) (low-bnd bnd))
		(setf (numeric-type-high res) (hi-bnd bnd))
		res))
	     (t
	      ;; Positive integer to a number is a number (for now)
	      (specifier-type 'number))))
	  ((csubtypep x-type (specifier-type 'rational))
	   ;; A rational to some power
1079
	   (case (numeric-type-class y-type)
rtoy's avatar
rtoy committed
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
	     (integer
	      ;; Positive rational to an integer power is always a rational
	      (specifier-type `(rational ,(or (low-bnd bnd) '*)
					 ,(or (hi-bnd bnd) '*))))
	     (rational
	      ;; Positive rational to rational power is either a rational
	      ;; or a single-float.
	      (let* ((lo (low-bnd bnd))
		     (hi (hi-bnd bnd))
		     (int-lo (if lo
				 (floor (bound-value lo))
				 '*))
		     (int-hi (if hi
				 (ceiling (bound-value hi))
				 '*))
		     (f-lo (if lo
			       (bound-func #'float lo)
			       '*))
		     (f-hi (if hi
			       (bound-func #'float hi)
			       '*)))
		(specifier-type `(or (rational ,int-lo ,int-hi)
				     (single-float ,f-lo, f-hi)))))
1103
	     (float
rtoy's avatar
rtoy committed
1104
1105
1106
1107
1108
	      ;; Positive rational to a float power is a float
	      (let ((res (copy-numeric-type y-type)))
		(setf (numeric-type-low res) (low-bnd bnd))
		(setf (numeric-type-high res) (hi-bnd bnd))
		res))
1109
	     (t
rtoy's avatar
rtoy committed
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
1137
1138
1139
1140
1141
1142
	      ;; Positive rational to a number is a number (for now)
	      (specifier-type 'number))))
	  ((csubtypep x-type (specifier-type 'float))
	   ;; A float to some power
	   (flet ((make-result (type)
		    (let ((res-type (or type 'float)))
		      (etypecase bnd
			(member-type
			 ;; Coerce all elements to the appropriate float
			 ;; type.
			 (make-member-type :members (mapcar #'(lambda (x)
								(coerce x res-type))
							    (member-type-members bnd))))
			(interval
			 (make-numeric-type
			  :class 'float
			  :format type
			  :low (coerce-numeric-bound (low-bnd bnd) res-type)
			  :high (coerce-numeric-bound (hi-bnd bnd) res-type)))))))
	     (case (numeric-type-class y-type)
	       ((or integer rational)
		;; Positive float to an integer or rational power is always a float
		(make-result (numeric-type-format x-type)))
	       (float
		;; Positive float to a float power is a float of the higher type
		(make-result (float-format-max (numeric-type-format x-type)
					       (numeric-type-format y-type))))
	       (t
		;; Positive float to a number is a number (for now)
		(specifier-type 'number)))))
	  (t
	   ;; A number to some power is a number.
	   (specifier-type 'number)))))
1143
  
1144
1145
1146
1147
1148
1149
1150
1151
1152
(defun merged-interval-expt (x y)
  (let* ((x-int (numeric-type->interval x))
	 (y-int (numeric-type->interval y)))
    (mapcar #'(lambda (type)
		(fixup-interval-expt type x-int y-int x y))
	    (flatten-list (interval-expt x-int y-int)))))

(defun expt-derive-type-aux (x y same-arg)
  (declare (ignore same-arg))
1153
1154
1155
  (cond ((or (not (numeric-type-real-p x))
	     (not (numeric-type-real-p y)))
	 ;; Use numeric contagion if either is not real
1156
1157
1158
1159
1160
	 (numeric-contagion x y))
	((csubtypep y (specifier-type 'integer))
	 ;; A real raised to an integer power is well-defined
	 (merged-interval-expt x y))
	(t
1161
	 ;; A real raised to a non-integral power is complicated....
1162
	 (cond ((or (csubtypep x (specifier-type '(rational 0)))
dtc's avatar
dtc committed
1163
		    (csubtypep x (specifier-type '(float (0d0)))))
1164
		;; A positive real to any power is well-defined.
1165
		(merged-interval-expt x y))
1166
1167
1168
1169
1170
	       ((and (csubtypep x (specifier-type 'rational))
		     (csubtypep x (specifier-type 'rational)))
		;; A rational to a rational power can be a rational or
		;; a single-float or a complex single-float.
		(specifier-type '(or rational single-float (complex single-float))))
1171
1172
1173
	       (t
		;; A real to some power.  The result could be a real
		;; or a complex.
1174
		(float-or-complex-float-type (numeric-contagion x y)))))))
1175
1176

(defoptimizer (expt derive-type) ((x y))
dtc's avatar
dtc committed
1177
  (two-arg-derive-type x y #'expt-derive-type-aux #'expt))
1178

ram's avatar
ram committed
1179

1180
;;; Note must assume that a type including 0.0 may also include -0.0
1181
;;; and thus the result may be complex -infinity + i*pi.
1182
;;;
1183
(defun log-derive-type-aux-1 (x)
1184
  (elfun-derive-type-simple x #'log 0d0 nil nil nil))
1185
1186

(defun log-derive-type-aux-2 (x y same-arg)
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
  (let ((log-x (log-derive-type-aux-1 x))
	(log-y (log-derive-type-aux-1 y))
	(result '()))
    ;; log-x or log-y might be union types.  We need to run through
    ;; the union types ourselves because /-derive-type-aux doesn't.
    (dolist (x-type (prepare-arg-for-derive-type log-x))
      (dolist (y-type (prepare-arg-for-derive-type log-y))
	(push (/-derive-type-aux x-type y-type same-arg) result)))
    (setf result (flatten-list result))
    (if (rest result)
	(make-union-type result)
	(first result))))
ram's avatar
ram committed
1199
1200

(defoptimizer (log derive-type) ((x &optional y))
1201
  (if y
dtc's avatar
dtc committed
1202
1203
      (two-arg-derive-type x y #'log-derive-type-aux-2 #'log)
      (one-arg-derive-type x #'log-derive-type-aux-1 #'log)))
1204
1205
1206


(defun atan-derive-type-aux-1 (y)
1207
  (elfun-derive-type-simple y #'atan nil nil (- (/ pi 2)) (/ pi 2)))
1208
1209
1210
1211

(defun atan-derive-type-aux-2 (y x same-arg)
  (declare (ignore same-arg))
  ;; The hard case with two args.  We just return the max bounds.
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
  (let ((result-type (numeric-contagion y x)))
    (cond ((and (numeric-type-real-p x)
		(numeric-type-real-p y))
	   (let* ((format (case (numeric-type-class result-type)
			    ((integer rational) 'single-float)
			    (t (numeric-type-format result-type))))
		  (bound-format (or format 'float)))
	     (make-numeric-type :class 'float
				:format format
				:complexp :real
				:low (coerce (- pi) bound-format)
				:high (coerce pi bound-format))))
	  (t
	   ;; The result is a float or a complex number
1226
	   (float-or-complex-float-type result-type)))))
ram's avatar
ram committed
1227
1228

(defoptimizer (atan derive-type) ((y &optional x))
1229
1230
1231
  (if x
      (two-arg-derive-type y x #'atan-derive-type-aux-2 #'atan)
      (one-arg-derive-type y #'atan-derive-type-aux-1 #'atan)))
1232
1233
1234


(defun cosh-derive-type-aux (x)
1235
  ;; We note that cosh x = cosh |x| for all real x.
1236
1237
1238
1239
  (elfun-derive-type-simple
   (if (numeric-type-real-p x)
       (abs-derive-type-aux x)
       x)
1240
   #'cosh nil nil 0 nil))
ram's avatar
ram committed
1241

1242
(defoptimizer (cosh derive-type) ((num))
dtc's avatar
dtc committed
1243
  (one-arg-derive-type num #'cosh-derive-type-aux #'cosh))
1244

1245

1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
(defun phase-derive-type-aux (arg)
  (let* ((format (case (numeric-type-class arg)
		   ((integer rational) 'single-float)
		   (t (numeric-type-format arg))))
	 (bound-type (or format 'float)))
    (cond ((numeric-type-real-p arg)
	   (case (interval-range-info (numeric-type->interval arg) 0.0)
	     ('+
	      ;; The number is positive, so the phase is 0.
	      (make-numeric-type :class 'float
				 :format format
				 :complexp :real
				 :low (coerce 0 bound-type)
				 :high (coerce 0 bound-type)))
	     ('-
	      ;; The number is always negative, so the phase is pi
	      (make-numeric-type :class 'float
				 :format format
				 :complexp :real
				 :low (coerce pi bound-type)
				 :high (coerce pi bound-type)))
	     (t
	      ;; We can't tell.  The result is 0 or pi.  Use a union
	      ;; type for this
	      (list
	       (make-numeric-type :class 'float
				  :format format
				  :complexp :real
				  :low (coerce 0 bound-type)
				  :high (coerce 0 bound-type))
	       (make-numeric-type :class 'float
				  :format format
				  :complexp :real
				  :low (coerce pi bound-type)
				  :high (coerce pi bound-type))))))
	  (t
	   ;; We have a complex number.  The answer is the range -pi
	   ;; to pi.  (-pi is included because we have -0.)
	   (make-numeric-type :class 'float
			      :format format
			      :complexp :real
			      :low (coerce (- pi) bound-type)
			      :high (coerce pi bound-type))))))
1289
1290

(defoptimizer (phase derive-type) ((num))
dtc's avatar
dtc committed
1291
  (one-arg-derive-type num #'phase-derive-type-aux #'phase))
dtc's avatar
dtc committed
1292

1293

1294
1295
1296
1297
(deftransform realpart ((x) ((complex rational)) *)
  '(kernel:%realpart x))
(deftransform imagpart ((x) ((complex rational)) *)
  '(kernel:%imagpart x))
ram's avatar
ram committed
1298

dtc's avatar
dtc committed
1299
1300
;;; Make REALPART and IMAGPART return the appropriate types.  This
;;; should help a lot in optimized code.
1301
1302

(defun realpart-derive-type-aux (type)
1303
1304
  (let ((class (numeric-type-class type))
	(format (numeric-type-format type)))
ram's avatar
ram committed
1305
    (cond ((numeric-type-real-p type)
1306
1307
	   ;; The realpart of a real has the same type and range as
	   ;; the input.
1308
1309
	   (make-numeric-type :class class
			      :format format
ram's avatar
ram committed
1310
1311
1312
1313
1314
			      :complexp :real
			      :low (numeric-type-low type)
			      :high (numeric-type-high type)))
	  (t
	   ;; We have a complex number.  The result has the same type
1315
	   ;; as the real part, except that it's real, not complex,
ram's avatar
ram committed
1316
	   ;; obviously.
1317
1318
	   (make-numeric-type :class class
			      :format format