/
float-tran.lisp
1414 lines (1289 loc) · 51.2 KB
/
float-tran.lisp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
;;;; This file contains floating-point-specific transforms, and may be
;;;; somewhat implementation-dependent in its assumptions of what the
;;;; formats are.
;;;; This software is part of the SBCL system. See the README file for
;;;; more information.
;;;;
;;;; This software is derived from the CMU CL system, which was
;;;; written at Carnegie Mellon University and released into the
;;;; public domain. The software is in the public domain and is
;;;; provided with absolutely no warranty. See the COPYING and CREDITS
;;;; files for more information.
(in-package "SB!C")
;;;; coercions
(defknown %single-float (real) single-float (movable foldable flushable))
(defknown %double-float (real) double-float (movable foldable flushable))
(deftransform float ((n f) (* single-float) *)
'(%single-float n))
(deftransform float ((n f) (* double-float) *)
'(%double-float n))
(deftransform float ((n) *)
'(if (floatp n)
n
(%single-float n)))
(deftransform %single-float ((n) (single-float) *)
'n)
(deftransform %double-float ((n) (double-float) *)
'n)
;;; RANDOM
(macrolet ((frob (fun type)
`(deftransform random ((num &optional state)
(,type &optional *) *)
"Use inline float operations."
'(,fun num (or state *random-state*)))))
(frob %random-single-float single-float)
(frob %random-double-float double-float))
;;; Mersenne Twister RNG
;;;
;;; FIXME: It's unpleasant to have RANDOM functionality scattered
;;; through the code this way. It would be nice to move this into the
;;; same file as the other RANDOM definitions.
(deftransform random ((num &optional state)
((integer 1 #.(expt 2 32)) &optional *))
;; FIXME: I almost conditionalized this as #!+sb-doc. Find some way
;; of automatically finding #!+sb-doc in proximity to DEFTRANSFORM
;; to let me scan for places that I made this mistake and didn't
;; catch myself.
"use inline (UNSIGNED-BYTE 32) operations"
(let ((num-high (numeric-type-high (lvar-type num))))
(when (null num-high)
(give-up-ir1-transform))
(cond ((constant-lvar-p num)
;; Check the worst case sum absolute 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 (- sb!kernel::random-integer-extra-bits)))
(give-up-ir1-transform
"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 (sb!bignum::%multiply
(random-chunk (or state *random-state*))
num)))))
((> num-high random-fixnum-max)
(give-up-ir1-transform
"The range is too large to ensure an accurate result."))
#!+x86
((< num-high (expt 2 32))
'(values (sb!bignum::%multiply (random-chunk (or state
*random-state*))
num)))
(t
'(rem (random-chunk (or state *random-state*)) num)))))
;;;; float accessors
(defknown make-single-float ((signed-byte 32)) single-float
(movable foldable flushable))
(defknown make-double-float ((signed-byte 32) (unsigned-byte 32)) double-float
(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))
(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)))
;;;; DECODE-FLOAT, INTEGER-DECODE-FLOAT, and SCALE-FLOAT
(defknown decode-single-float (single-float)
(values single-float single-float-exponent (single-float -1f0 1f0))
(movable foldable flushable))
(defknown decode-double-float (double-float)
(values double-float double-float-exponent (double-float -1d0 1d0))
(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))
(movable foldable flushable))
(defknown scale-single-float (single-float fixnum) single-float
(movable foldable flushable))
(defknown scale-double-float (double-float fixnum) double-float
(movable foldable flushable))
(deftransform decode-float ((x) (single-float) *)
'(decode-single-float x))
(deftransform decode-float ((x) (double-float) *)
'(decode-double-float x))
(deftransform integer-decode-float ((x) (single-float) *)
'(integer-decode-single-float x))
(deftransform integer-decode-float ((x) (double-float) *)
'(integer-decode-double-float x))
(deftransform scale-float ((f ex) (single-float *) *)
(if (and #!+x86 t #!-x86 nil
(csubtypep (lvar-type ex)
(specifier-type '(signed-byte 32))))
'(coerce (%scalbn (coerce f 'double-float) ex) 'single-float)
'(scale-single-float f ex)))
(deftransform scale-float ((f ex) (double-float *) *)
(if (and #!+x86 t #!-x86 nil
(csubtypep (lvar-type ex)
(specifier-type '(signed-byte 32))))
'(%scalbn f ex)
'(scale-double-float f ex)))
;;; What is the CROSS-FLOAT-INFINITY-KLUDGE?
;;;
;;; SBCL's own implementation of floating point supports floating
;;; point infinities. Some of the old CMU CL :PROPAGATE-FLOAT-TYPE and
;;; :PROPAGATE-FUN-TYPE code, like the DEFOPTIMIZERs below, uses this
;;; floating point support. Thus, we have to avoid running it on the
;;; cross-compilation host, since we're not guaranteed that the
;;; cross-compilation host will support floating point infinities.
;;;
;;; If we wanted to live dangerously, we could conditionalize the code
;;; with #+(OR SBCL SB-XC) instead. That way, if the cross-compilation
;;; host happened to be SBCL, we'd be able to run the infinity-using
;;; code. Pro:
;;; * SBCL itself gets built with more complete optimization.
;;; Con:
;;; * You get a different SBCL depending on what your cross-compilation
;;; host is.
;;; So far the pros and cons seem seem to be mostly academic, since
;;; AFAIK (WHN 2001-08-28) the propagate-foo-type optimizations aren't
;;; actually important in compiling SBCL itself. If this changes, then
;;; we have to decide:
;;; * Go for simplicity, leaving things as they are.
;;; * Go for performance at the expense of conceptual clarity,
;;; using #+(OR SBCL SB-XC) and otherwise leaving the build
;;; process as is.
;;; * Go for performance at the expense of build time, using
;;; #+(OR SBCL SB-XC) and also making SBCL do not just
;;; make-host-1.sh and make-host-2.sh, but a third step
;;; make-host-3.sh where it builds itself under itself. (Such a
;;; 3-step build process could also help with other things, e.g.
;;; using specialized arrays to represent debug information.)
;;; * Rewrite the code so that it doesn't depend on unportable
;;; floating point infinities.
;;; optimizers for SCALE-FLOAT. If the float has bounds, new bounds
;;; are computed for the result, if possible.
#-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
(progn
(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
(handler-case
(scale-float (type-bound-number x) n)
(floating-point-overflow ()
nil))
(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)))
(make-numeric-type :class (numeric-type-class f)
:format (numeric-type-format f)
:complexp :real
:low new-lo
:high new-hi)))))
(defoptimizer (scale-single-float derive-type) ((f ex))
(two-arg-derive-type f ex #'scale-float-derive-type-aux
#'scale-single-float t))
(defoptimizer (scale-double-float derive-type) ((f ex))
(two-arg-derive-type f ex #'scale-float-derive-type-aux
#'scale-double-float t))
;;; DEFOPTIMIZERs for %SINGLE-FLOAT and %DOUBLE-FLOAT. This makes the
;;; FLOAT function return the correct ranges if the input has some
;;; defined range. Quite useful if we want to convert some type of
;;; bounded integer into a float.
(macrolet
((frob (fun type)
(let ((aux-name (symbolicate fun "-DERIVE-TYPE-AUX")))
`(progn
(defun ,aux-name (num)
;; When converting a number to a float, the limits are
;; the same.
(let* ((lo (bound-func (lambda (x)
(coerce x ',type))
(numeric-type-low num)))
(hi (bound-func (lambda (x)
(coerce x ',type))
(numeric-type-high num))))
(specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
(defoptimizer (,fun derive-type) ((num))
(one-arg-derive-type num #',aux-name #',fun))))))
(frob %single-float single-float)
(frob %double-float double-float))
) ; PROGN
;;;; float contagion
;;; Do some stuff to recognize when the loser is doing mixed float and
;;; rational arithmetic, or different float types, and fix it up. If
;;; we don't, he won't even get so much as an efficiency note.
(deftransform float-contagion-arg1 ((x y) * * :defun-only t :node node)
`(,(lvar-fun-name (basic-combination-fun node))
(float x y) y))
(deftransform float-contagion-arg2 ((x y) * * :defun-only t :node node)
`(,(lvar-fun-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))
(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, and 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).
(macrolet ((frob (op)
`(deftransform ,op ((x y) (float rational) *)
"open-code FLOAT to RATIONAL comparison"
(unless (constant-lvar-p y)
(give-up-ir1-transform
"The RATIONAL value isn't known at compile time."))
(let ((val (lvar-value y)))
(unless (eql (rational (float val)) val)
(give-up-ir1-transform
"~S doesn't have a precise float representation."
val)))
`(,',op x (float y x)))))
(frob <)
(frob >)
(frob =))
;;;; irrational derive-type methods
;;; Derive the result to be float for argument types in the
;;; appropriate domain.
#+sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
(dolist (stuff '((asin (real -1.0 1.0))
(acos (real -1.0 1.0))
(acosh (real 1.0))
(atanh (real -1.0 1.0))
(sqrt (real 0.0))))
(destructuring-bind (name type) stuff
(let ((type (specifier-type type)))
(setf (fun-info-derive-type (fun-info-or-lose name))
(lambda (call)
(declare (type combination call))
(when (csubtypep (lvar-type
(first (combination-args call)))
type)
(specifier-type 'float)))))))
#+sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
(defoptimizer (log derive-type) ((x &optional y))
(when (and (csubtypep (lvar-type x)
(specifier-type '(real 0.0)))
(or (null y)
(csubtypep (lvar-type y)
(specifier-type '(real 0.0)))))
(specifier-type 'float)))
;;;; irrational transforms
(defknown (%tan %sinh %asinh %atanh %log %logb %log10 %tan-quick)
(double-float) double-float
(movable foldable flushable))
(defknown (%sin %cos %tanh %sin-quick %cos-quick)
(double-float) (double-float -1.0d0 1.0d0)
(movable foldable flushable))
(defknown (%asin %atan)
(double-float)
(double-float #.(coerce (- (/ pi 2)) 'double-float)
#.(coerce (/ pi 2) 'double-float))
(movable foldable flushable))
(defknown (%acos)
(double-float) (double-float 0.0d0 #.(coerce pi 'double-float))
(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))
(defknown %expm1
(double-float) (double-float -1d0)
(movable foldable flushable))
(defknown (%hypot)
(double-float double-float) (double-float 0d0)
(movable foldable flushable))
(defknown (%pow)
(double-float double-float) double-float
(movable foldable flushable))
(defknown (%atan2)
(double-float double-float)
(double-float #.(coerce (- pi) 'double-float)
#.(coerce pi 'double-float))
(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))
(defknown (%log1p)
(double-float) double-float
(movable foldable flushable))
(macrolet ((def (name prim rtype)
`(progn
(deftransform ,name ((x) (single-float) ,rtype)
`(coerce (,',prim (coerce x 'double-float)) 'single-float))
(deftransform ,name ((x) (double-float) ,rtype)
`(,',prim x)))))
(def exp %exp *)
(def log %log float)
(def sqrt %sqrt float)
(def asin %asin float)
(def acos %acos float)
(def atan %atan *)
(def sinh %sinh *)
(def cosh %cosh *)
(def tanh %tanh *)
(def asinh %asinh *)
(def acosh %acosh float)
(def atanh %atanh float))
;;; 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.
(macrolet ((def (name prim prim-quick)
(declare (ignorable prim-quick))
`(progn
(deftransform ,name ((x) (single-float) *)
#!+x86 (cond ((csubtypep (lvar-type x)
(specifier-type '(single-float
(#.(- (expt 2f0 64)))
(#.(expt 2f0 64)))))
`(coerce (,',prim-quick (coerce x 'double-float))
'single-float))
(t
(compiler-notify
"unable to avoid inline argument range check~@
because the argument range (~S) was not within 2^64"
(type-specifier (lvar-type x)))
`(coerce (,',prim (coerce x 'double-float)) 'single-float)))
#!-x86 `(coerce (,',prim (coerce x 'double-float)) 'single-float))
(deftransform ,name ((x) (double-float) *)
#!+x86 (cond ((csubtypep (lvar-type x)
(specifier-type '(double-float
(#.(- (expt 2d0 64)))
(#.(expt 2d0 64)))))
`(,',prim-quick x))
(t
(compiler-notify
"unable to avoid inline argument range check~@
because the argument range (~S) was not within 2^64"
(type-specifier (lvar-type x)))
`(,',prim x)))
#!-x86 `(,',prim x)))))
(def sin %sin %sin-quick)
(def cos %cos %cos-quick)
(def tan %tan %tan-quick))
(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) *)
`(%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) *)
`(%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)) *)
`(%pow x (coerce y 'double-float)))
;;; ANSI says log with base zero returns zero.
(deftransform log ((x y) (float float) float)
'(if (zerop y) y (/ (log x) (log y))))
;;; Handle some simple transformations.
(deftransform abs ((x) ((complex double-float)) double-float)
'(%hypot (realpart x) (imagpart x)))
(deftransform abs ((x) ((complex single-float)) single-float)
'(coerce (%hypot (coerce (realpart x) 'double-float)
(coerce (imagpart x) 'double-float))
'single-float))
(deftransform phase ((x) ((complex double-float)) double-float)
'(%atan2 (imagpart x) (realpart x)))
(deftransform phase ((x) ((complex single-float)) single-float)
'(coerce (%atan2 (coerce (imagpart x) 'double-float)
(coerce (realpart x) 'double-float))
'single-float))
(deftransform phase ((x) ((float)) float)
'(if (minusp (float-sign x))
(float pi x)
(float 0 x)))
;;; The number is of type REAL.
(defun numeric-type-real-p (type)
(and (numeric-type-p type)
(eq (numeric-type-complexp type) :real)))
;;; Coerce a numeric type bound to the given type while handling
;;; exclusive bounds.
(defun coerce-numeric-bound (bound type)
(when bound
(if (consp bound)
(list (coerce (car bound) type))
(coerce bound type))))
#-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
(progn
;;;; optimizers for elementary functions
;;;;
;;;; These optimizers compute the output range of the elementary
;;;; function, based on the domain of the input.
;;; 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))))
;;; Compute a specifier like '(OR FLOAT (COMPLEX FLOAT)), except float
;;; should be the right kind of float. Allow bounds for the float
;;; part too.
(defun float-or-complex-float-type (arg &optional lo hi)
(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))
(lo (coerce-numeric-bound lo float-type))
(hi (coerce-numeric-bound hi float-type)))
(specifier-type `(or (,float-type ,(or lo '*) ,(or hi '*))
(complex ,float-type)))))
) ; PROGN
(eval-when (:compile-toplevel :execute)
;; So the problem with this hack is that it's actually broken. If
;; the host does not have long floats, then setting *R-D-F-F* to
;; LONG-FLOAT doesn't actually buy us anything. FIXME.
(setf *read-default-float-format*
#!+long-float 'long-float #!-long-float 'double-float))
;;; Test whether the numeric-type ARG is within in domain specified by
;;; DOMAIN-LOW and DOMAIN-HIGH, consider negative and positive zero to
;;; be distinct.
#-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
(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 (type-bound-number arg-lo))
(arg-hi (numeric-type-high arg))
(arg-hi-val (type-bound-number arg-hi)))
;; Check that the ARG bounds are correctly canonicalized.
(when (and arg-lo (floatp arg-lo-val) (zerop arg-lo-val) (consp arg-lo)
(minusp (float-sign arg-lo-val)))
(compiler-notify "float zero bound ~S not correctly canonicalized?" arg-lo)
(setq arg-lo 0e0 arg-lo-val arg-lo))
(when (and arg-hi (zerop arg-hi-val) (floatp arg-hi-val) (consp arg-hi)
(plusp (float-sign arg-hi-val)))
(compiler-notify "float zero bound ~S not correctly canonicalized?" arg-hi)
(setq arg-hi (ecase *read-default-float-format*
(double-float (load-time-value (make-unportable-float :double-float-negative-zero)))
#!+long-float
(long-float (load-time-value (make-unportable-float :long-float-negative-zero))))
arg-hi-val arg-hi))
(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)))))))))
(eval-when (:compile-toplevel :execute)
(setf *read-default-float-format* 'single-float))
#-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
(progn
;;; Handle monotonic functions of a single variable whose domain is
;;; possibly part of the real line. ARG is the variable, FCN is the
;;; 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.
;;;
;;; Negative and positive zero are considered distinct within
;;; DOMAIN-LOW and DOMAIN-HIGH.
;;;
;;; DEFAULT-LOW and DEFAULT-HIGH are the lower and upper bounds if we
;;; can't compute the bounds using FCN.
(defun elfun-derive-type-simple (arg fcn domain-low domain-high
default-low default-high
&optional (increasingp t))
(declare (type (or null real) domain-low domain-high))
(etypecase arg
(numeric-type
(cond ((eq (numeric-type-complexp arg) :complex)
(make-numeric-type :class (numeric-type-class arg)
:format (numeric-type-format arg)
:complexp :complex))
((numeric-type-real-p arg)
;; 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)
(make-interval
:low domain-low
:high domain-high))
(cond
(intersection
;; Process the intersection.
(let* ((low (interval-low intersection))
(high (interval-high intersection))
(res-lo (or (bound-func fcn (if increasingp low high))
default-low))
(res-hi (or (bound-func fcn (if increasingp high low))
default-high))
(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.
(if (or (null difference)
;; Check whether the arg is within the domain.
(domain-subtypep arg domain-low domain-high))
result-type
(list result-type
(specifier-type `(complex ,bound-type))))))
(t
;; No intersection so the result must be purely complex.
(complex-float-type arg)))))
(t
(float-or-complex-float-type arg default-low default-high))))))
(macrolet
((frob (name domain-low domain-high def-low-bnd def-high-bnd
&key (increasingp t))
(let ((num (gensym)))
`(defoptimizer (,name derive-type) ((,num))
(one-arg-derive-type
,num
(lambda (arg)
(elfun-derive-type-simple arg #',name
,domain-low ,domain-high
,def-low-bnd ,def-high-bnd
,increasingp))
#',name)))))
;; These functions are easy because they are defined for the whole
;; real line.
(frob exp nil nil 0 nil)
(frob sinh nil nil nil nil)
(frob tanh nil nil -1 1)
(frob asinh nil nil nil nil)
;; These functions are only defined for part of the real line. The
;; condition selects the desired part of the line.
(frob asin -1d0 1d0 (- (/ pi 2)) (/ pi 2))
;; Acos is monotonic decreasing, so we need to swap the function
;; values at the lower and upper bounds of the input domain.
(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 (load-time-value (make-unportable-float :double-float-negative-zero)) nil 0 nil))
;;; 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.
(defun safe-expt (x y)
(handler-case
(when (< (abs y) 10000)
(expt x y))
(error ()
nil)))
;;; Handle the case when x >= 1.
(defun interval-expt-> (x y)
(case (sb!c::interval-range-info y 0d0)
(+
;; Y is positive and log X >= 0. The range of exp(y * log(x)) is
;; obviously non-negative. We just have to be careful for
;; infinite bounds (given by nil).
(let ((lo (safe-expt (type-bound-number (sb!c::interval-low x))
(type-bound-number (sb!c::interval-low y))))
(hi (safe-expt (type-bound-number (sb!c::interval-high x))
(type-bound-number (sb!c::interval-high y)))))
(list (sb!c::make-interval :low (or lo 1) :high hi))))
(-
;; Y is negative and log x >= 0. The range of exp(y * log(x)) is
;; obviously [0, 1]. However, underflow (nil) means 0 is the
;; result.
(let ((lo (safe-expt (type-bound-number (sb!c::interval-high x))
(type-bound-number (sb!c::interval-low y))))
(hi (safe-expt (type-bound-number (sb!c::interval-low x))
(type-bound-number (sb!c::interval-high y)))))
(list (sb!c::make-interval :low (or lo 0) :high (or hi 1)))))
(t
;; Split the interval in half.
(destructuring-bind (y- y+)
(sb!c::interval-split 0 y t)
(list (interval-expt-> x y-)
(interval-expt-> x y+))))))
;;; Handle the case when x <= 1
(defun interval-expt-< (x y)
(case (sb!c::interval-range-info x 0d0)
(+
;; The case of 0 <= x <= 1 is easy
(case (sb!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 (type-bound-number (sb!c::interval-low x))
(type-bound-number (sb!c::interval-high y))))
(hi (safe-expt (type-bound-number (sb!c::interval-high x))
(type-bound-number (sb!c::interval-low y)))))
(list (sb!c::make-interval :low (or lo 0) :high (or hi 1)))))
(-
;; Y is negative and log x <= 0. The range of exp(y * log(x)) is
;; obviously [1, inf].
(let ((hi (safe-expt (type-bound-number (sb!c::interval-low x))
(type-bound-number (sb!c::interval-low y))))
(lo (safe-expt (type-bound-number (sb!c::interval-high x))
(type-bound-number (sb!c::interval-high y)))))
(list (sb!c::make-interval :low (or lo 1) :high hi))))
(t
;; Split the interval in half
(destructuring-bind (y- y+)
(sb!c::interval-split 0 y t)
(list (interval-expt-< x y-)
(interval-expt-< x y+))))))
(-
;; 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 (sb!c::make-interval :low nil :high nil)))
(t
(destructuring-bind (neg pos)
(interval-split 0 x t t)
(list (interval-expt-< neg y)
(interval-expt-< pos y))))))
;;; Compute bounds for (expt x y).
(defun interval-expt (x y)
(case (interval-range-info x 1)
(+
;; X >= 1
(interval-expt-> x y))
(-
;; X <= 1
(interval-expt-< x y))
(t
(destructuring-bind (left right)
(interval-split 1 x t t)
(list (interval-expt left y)
(interval-expt right y))))))
(defun fixup-interval-expt (bnd x-int y-int x-type y-type)
(declare (ignore x-int))
;; Figure out what the return type should be, given the argument
;; types and bounds and the result type and bounds.
(cond ((csubtypep x-type (specifier-type 'integer))
;; an integer to some power
(case (numeric-type-class y-type)
(integer
;; Positive integer to an integer power is either an
;; integer or a rational.
(let ((lo (or (interval-low bnd) '*))
(hi (or (interval-high bnd) '*)))
(if (and (interval-low y-int)
(>= (type-bound-number (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 (interval-low bnd))
(hi (interval-high bnd))
(int-lo (if lo
(floor (type-bound-number lo))
'*))
(int-hi (if hi
(ceiling (type-bound-number 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)))))
(float
;; A positive integer to a float power is a float.
(modified-numeric-type y-type
:low (interval-low bnd)
:high (interval-high bnd)))
(t
;; A positive integer to a number is a number (for now).
(specifier-type 'number))))
((csubtypep x-type (specifier-type 'rational))
;; a rational to some power
(case (numeric-type-class y-type)
(integer
;; A positive rational to an integer power is always a rational.
(specifier-type `(rational ,(or (interval-low bnd) '*)
,(or (interval-high bnd) '*))))
(rational
;; A positive rational to rational power is either a rational
;; or a single-float.
(let* ((lo (interval-low bnd))
(hi (interval-high bnd))
(int-lo (if lo
(floor (type-bound-number lo))
'*))
(int-hi (if hi
(ceiling (type-bound-number 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)))))
(float
;; A positive rational to a float power is a float.
(modified-numeric-type y-type
:low (interval-low bnd)
:high (interval-high bnd)))
(t
;; A positive rational to a number is a number (for now).
(specifier-type 'number))))
((csubtypep x-type (specifier-type 'float))
;; a float to some power
(case (numeric-type-class y-type)
((or integer rational)
;; A positive float to an integer or rational power is
;; always a float.
(make-numeric-type
:class 'float
:format (numeric-type-format x-type)
:low (interval-low bnd)
:high (interval-high bnd)))
(float
;; A positive float to a float power is a float of the
;; higher type.
(make-numeric-type
:class 'float
:format (float-format-max (numeric-type-format x-type)
(numeric-type-format y-type))
:low (interval-low bnd)
:high (interval-high bnd)))
(t
;; A 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))))
(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))
(cond ((or (not (numeric-type-real-p x))
(not (numeric-type-real-p y)))
;; Use numeric contagion if either is not real.
(numeric-contagion x y))
((csubtypep y (specifier-type 'integer))
;; A real raised to an integer power is well-defined.
(merged-interval-expt x y))
;; A real raised to a non-integral power can be a float or a
;; complex number.
((or (csubtypep x (specifier-type '(rational 0)))
(csubtypep x (specifier-type '(float (0d0)))))
;; But a positive real to any power is well-defined.
(merged-interval-expt x y))
((and (csubtypep x (specifier-type 'rational))
(csubtypep x (specifier-type 'rational)))
;; A rational to the power of a rational could be a rational
;; or a possibly-complex single float
(specifier-type '(or rational single-float (complex single-float))))
(t
;; a real to some power. The result could be a real or a
;; complex.
(float-or-complex-float-type (numeric-contagion x y)))))
(defoptimizer (expt derive-type) ((x y))
(two-arg-derive-type x y #'expt-derive-type-aux #'expt))
;;; Note we must assume that a type including 0.0 may also include
;;; -0.0 and thus the result may be complex -infinity + i*pi.
(defun log-derive-type-aux-1 (x)
(elfun-derive-type-simple x #'log 0d0 nil nil nil))
(defun log-derive-type-aux-2 (x y same-arg)
(let ((log-x (log-derive-type-aux-1 x))
(log-y (log-derive-type-aux-1 y))
(accumulated-list nil))
;; 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) accumulated-list)))
(apply #'type-union (flatten-list accumulated-list))))
(defoptimizer (log derive-type) ((x &optional y))
(if y
(two-arg-derive-type x y #'log-derive-type-aux-2 #'log)
(one-arg-derive-type x #'log-derive-type-aux-1 #'log)))
(defun atan-derive-type-aux-1 (y)
(elfun-derive-type-simple y #'atan nil nil (- (/ pi 2)) (/ pi 2)))
(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.
(let ((result-type (numeric-contagion y x)))
(cond ((and (numeric-type-real-p x)
(numeric-type-real-p y))
(let* (;; FIXME: This expression for FORMAT seems to
;; appear multiple times, and should be factored out.
(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
(float-or-complex-float-type result-type)))))
(defoptimizer (atan derive-type) ((y &optional x))
(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)))
(defun cosh-derive-type-aux (x)
;; We note that cosh x = cosh |x| for all real x.
(elfun-derive-type-simple
(if (numeric-type-real-p x)
(abs-derive-type-aux x)
x)
#'cosh nil nil 0 nil))
(defoptimizer (cosh derive-type) ((num))
(one-arg-derive-type num #'cosh-derive-type-aux #'cosh))
(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