source: trunk/source/level-0/l0-float.lisp @ 14365

Last change on this file since 14365 was 14365, checked in by gb, 9 years ago

Admit that the macro CCL::REPORT-BAD-ARG has accepted exactly 2 arguments
for the last 20+ years. If the second argument (the typespec) is quoted,
warn at macroexpand time if it looks suspicious.

Fix a handful of cases that were detected by that change. In one
case, introduce a predicate so that EXPT can complain about an
argument whose realpart isn't positive via a SATISFIES type specifier.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.3 KB
Line 
1;;; -*- Mode: Lisp; Package: CCL -*-
2;;;
3;;;   Copyright (C) 2009 Clozure Associates
4;;;   Copyright (C) 1994-2001 Digitool, Inc
5;;;   This file is part of Clozure CL. 
6;;;
7;;;   Clozure CL is licensed under the terms of the Lisp Lesser GNU Public
8;;;   License , known as the LLGPL and distributed with Clozure CL as the
9;;;   file "LICENSE".  The LLGPL consists of a preamble and the LGPL,
10;;;   which is distributed with Clozure CL as the file "LGPL".  Where these
11;;;   conflict, the preamble takes precedence. 
12;;;
13;;;   Clozure CL is referenced in the preamble as the "LIBRARY."
14;;;
15;;;   The LLGPL is also available online at
16;;;   http://opensource.franz.com/preamble.html
17;;;
18;;; level-0;l0-float.lisp
19
20(in-package "CCL")
21
22(eval-when (:compile-toplevel :execute)
23  (require "NUMBER-MACROS")
24  (require :number-case-macro) 
25)
26
27;;; used by float reader
28(defun make-float-from-fixnums (hi lo exp sign &optional result)
29  ;;(require-null-or-double-float-sym result)
30  ;; maybe nuke all these require-types?
31  ;;(setq hi (require-type hi 'fixnum))
32  ;;(setq lo (require-type lo 'fixnum))
33  ;;(setq exp (require-type exp 'fixnum))
34  ;;(setq sign (require-type sign 'fixnum))
35  (let ((the-float (or result (%make-dfloat))))
36    (%make-float-from-fixnums the-float hi lo exp sign)
37    the-float))
38
39
40#+32-bit-target
41(defun make-short-float-from-fixnums (significand biased-exp sign &optional result)
42  (%make-short-float-from-fixnums (or result (%make-sfloat)) significand biased-exp sign))
43
44#+64-bit-target
45(defun make-short-float-from-fixnums (significand biased-exp sign)
46  (declare (fixnum significand biased-exp sign))
47  (host-single-float-from-unsigned-byte-32
48   (logior
49    (the fixnum (if (< sign 0) (ash 1 31) 0))
50    (the fixnum (ash biased-exp IEEE-single-float-exponent-offset))
51    (the fixnum (logand significand
52                        (1- (ash 1 IEEE-single-float-hidden-bit)))))))
53
54
55(defun float-sign (n1 &optional n2) ; second arg silly
56  "Return a floating-point number that has the same sign as
57   FLOAT1 and, if FLOAT2 is given, has the same absolute value
58   as FLOAT2."
59  (if (and n2 (not (typep n2 'float)))
60    (setq n2 (require-type n2 'float)))
61  (number-case n1
62    (double-float                       
63     (if (%double-float-sign n1) 
64       (if n2
65         (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) n2 (- n2))
66         -1.0d0)
67       (if n2
68         (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) (- n2) n2)
69         1.0d0)))
70    (short-float
71     (if (%short-float-sign n1)
72       (if n2
73         (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) n2 (- n2))
74         -1.0s0)
75       (if n2
76         (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) (- n2) n2)
77         1.0s0)))))
78
79
80
81(defun %double-float-minusp (n)
82  (and (%double-float-sign n)(not (%double-float-zerop n))))
83
84(defun %short-float-minusp (n)
85  (and (%short-float-sign n) (not (%short-float-zerop n))))
86
87(defun %double-float-abs (n)
88  (if (not (%double-float-sign n))
89    n 
90    (%%double-float-abs! n (%make-dfloat))))
91
92#+32-bit-target
93(defun %short-float-abs (n)
94  (if (not (%short-float-sign n))
95    n 
96    (%%short-float-abs! n (%make-sfloat))))
97
98(defun fixnum-decode-float (n)
99  (etypecase n
100    (double-float (%integer-decode-double-float n))))
101
102(defun nan-or-infinity-p (n)
103  (etypecase n
104    (double-float (eq 2047 (%double-float-exp n)))
105    (short-float (eq 255 (%short-float-exp n)))))
106
107; not sure this is right
108(defun infinity-p (n)
109  (etypecase n
110    (double-float (multiple-value-bind (hi lo exp)(fixnum-decode-float n)
111                    (and (eq 2047 exp)
112                         (eq #x1000000 hi)
113                         (eq 0 lo))))
114    (short-float
115     #+32-bit-target
116     (multiple-value-bind (high low)(%sfloat-hwords n)
117                  (let*  ((mantissa (%ilogior2 low (%ilsl 16 (%ilogand2 high #x007F))))
118                          (exp (%ilsr 7 (%ilogand2 high #x7F80))))
119                    (and (eq exp 255)
120                         (eq 0 mantissa))))
121     #+64-bit-target
122     (let* ((bits (single-float-bits n))
123            (exp (ldb (byte IEEE-single-float-exponent-width
124                            IEEE-single-float-exponent-offset)
125                      bits))
126            (mantissa (ldb (byte IEEE-single-float-mantissa-width
127                            IEEE-single-float-mantissa-offset)
128                           bits)))
129       (declare (fixnum bits exp mantissa))
130       (and (= exp 255)
131            (zerop mantissa))))))
132
133#+32-bit-target
134(defun fixnum-decode-short-float (float)
135  (multiple-value-bind (high low)(%sfloat-hwords float)
136    (let*  ((mantissa (%ilogior2 low (%ilsl 16 (%ilogand2 high #x007F))))
137            (exp (%ilsr 7 (%ilogand2 high #x7F80))))
138      (if (and (neq exp 0)(neq exp 255))(setq mantissa (%ilogior mantissa #x800000)))
139      (values mantissa exp (%ilsr 15 high)))))
140
141#+64-bit-target
142(defun fixnum-decode-short-float (float)
143  (let* ((bits (single-float-bits float)))
144    (declare (fixnum bits))
145    (let* ((mantissa (ldb (byte IEEE-single-float-mantissa-width
146                                IEEE-single-float-mantissa-offset)
147                          bits))
148           (exp (ldb (byte IEEE-single-float-exponent-width
149                           IEEE-single-float-exponent-offset)
150                     bits))
151           (sign (lsh bits -31)))
152      (declare (fixnum mantissa exp sign))
153      (unless (or (= exp 0) (= exp 255))
154        (setq mantissa (logior mantissa (ash 1 IEEE-single-float-hidden-bit))))
155      (values mantissa exp sign))))
156                 
157                   
158
159#+32-bit-target
160(defun integer-decode-double-float (n)
161  (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n)
162    ; is only 53 bits and positive so should be easy
163    ;(values (logior (ash hi 28) lo) exp sign)))
164    ; if denormalized, may fit in a fixnum
165    (setq exp (- exp (if (< hi #x1000000) 
166                       (+ IEEE-double-float-mantissa-width IEEE-double-float-bias)
167                       (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias)))))
168    (if (< hi (ash 1 (1- target::fixnumshift))) ; aka 2
169      (values (logior (ash hi 28) lo) exp sign)
170      ; might fit in 1 word?
171      (let ((big (%alloc-misc 2 target::subtag-bignum)))
172        (make-big-53 hi lo big)
173        (if (< hi #x1000000) (%normalize-bignum big))
174        (values big exp sign)))))
175
176#+64-bit-target
177(defun integer-decode-double-float (n)
178  (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n)
179    (setq exp (- exp (if (< hi #x1000000) 
180                       (+ IEEE-double-float-mantissa-width IEEE-double-float-bias)
181                       (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias)))))
182    (values (logior (ash hi 28) lo) exp sign)))
183   
184
185;;; actually only called when magnitude bigger than a fixnum
186#+32-bit-target
187(defun %truncate-double-float (n)
188  (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n)
189    (if (< exp (1+ IEEE-double-float-bias)) ; this is false in practice
190      0
191      (progn
192        (setq exp (- exp (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias))))
193        (if (eq sign 1)  ; positive
194          (logior (ash hi (+ 28 exp))(ash lo exp))
195          (if (<= exp 0) ; exp positive - negate before shift - else after
196            (let ((poo (logior (ash hi (+ 28 exp))(ash lo exp))))
197              (- poo))
198            (let ((poo (logior (ash hi 28) lo)))
199              (ash (- poo) exp))))))))
200
201#+64-bit-target
202(defun %truncate-double-float (n)
203  (multiple-value-bind (mantissa exp sign) (integer-decode-float n)
204    (* sign (ash mantissa exp))))
205
206
207
208; actually only called when bigger than a fixnum
209(defun %truncate-short-float (n)
210  (multiple-value-bind (mantissa exp sign)(fixnum-decode-short-float n)
211    (if (< exp (1+ IEEE-single-float-bias)) ; is magnitude less than 1 - false in practice
212      0
213      (progn
214        (setq exp (- exp (+ IEEE-single-float-mantissa-width (1+ IEEE-single-float-bias))))
215        (ash (if (eq sign 0) mantissa (- mantissa)) exp)))))
216
217(defun decode-float (n)
218  "Return three values:
219   1) a floating-point number representing the significand. This is always
220      between 0.5 (inclusive) and 1.0 (exclusive).
221   2) an integer representing the exponent.
222   3) -1.0 or 1.0 (i.e. the sign of the argument.)"
223  (number-case n
224    (double-float
225     (let* ((old-exp (%double-float-exp n))
226            (sign (if (%double-float-sign n) -1.0d0 1.0d0)))   
227       (if (eq 0 old-exp)
228         (if (%double-float-zerop n)
229           (values 0.0d0 0 sign)
230           (let* ((val (%make-dfloat))
231                  (zeros (dfloat-significand-zeros n)))
232             (%%double-float-abs! n val)
233             (%%scale-dfloat! val (+ 2 IEEE-double-float-bias zeros) val) ; get it normalized
234             (set-%double-float-exp val IEEE-double-float-bias)      ; then bash exponent
235             (values val (- old-exp zeros IEEE-double-float-bias) sign)))
236         (if (> old-exp IEEE-double-float-normal-exponent-max)
237           (error "Can't decode NAN or infinity ~s" n)
238           (let ((val (%make-dfloat)))
239             (%%double-float-abs! n val)
240             (set-%double-float-exp val IEEE-double-float-bias)
241             (values val (- old-exp IEEE-double-float-bias) sign))))))
242    (short-float
243     (let* ((old-exp (%short-float-exp n))
244            (sign (if (%short-float-sign n) -1.0s0 1.0s0)))
245       (if (eq 0 old-exp)
246         (if (%short-float-zerop n)
247           (values 0.0s0 0 sign)
248           #+32-bit-target
249           (let* ((val (%make-sfloat))
250                  (zeros (sfloat-significand-zeros n)))
251             (%%short-float-abs! n val)
252             (%%scale-sfloat! val (+ 2 IEEE-single-float-bias zeros) val) ; get it normalized
253             (set-%short-float-exp val IEEE-single-float-bias)      ; then bash exponent
254             (values val (- old-exp zeros IEEE-single-float-bias) sign))
255           #+64-bit-target
256           (let* ((zeros (sfloat-significand-zeros n))
257                  (val (%%scale-sfloat (%short-float-abs n)
258                                       (+ 2 IEEE-single-float-bias zeros))))
259             (values (set-%short-float-exp val IEEE-single-float-bias)
260                     (- old-exp zeros IEEE-single-float-bias) sign)))
261         (if (> old-exp IEEE-single-float-normal-exponent-max)
262           (error "Can't decode NAN or infinity ~s" n)
263           #+32-bit-target
264           (let ((val (%make-sfloat)))
265             (%%short-float-abs! n val)
266             (set-%short-float-exp val IEEE-single-float-bias)
267             (values val (- old-exp IEEE-single-float-bias) sign))
268           #+64-bit-target
269           (values (set-%short-float-exp (%short-float-abs n)
270                                         IEEE-single-float-bias)
271                   (- old-exp IEEE-single-float-bias) sign)))))))
272
273; (* float (expt 2 int))
274(defun scale-float (float int)
275  "Return the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
276  of precision or overflow."
277  (unless (fixnump int)(setq int (require-type int 'fixnum)))
278  (number-case float
279    (double-float
280     (let* ((float-exp (%double-float-exp float))
281            (new-exp (+ float-exp int)))
282       (if (eq 0 float-exp) ; already denormalized?
283         (if (%double-float-zerop float)
284           float 
285           (let ((result (%make-dfloat)))
286             (%%scale-dfloat! float (+ (1+ IEEE-double-float-bias) int) result)))
287         (if (<= new-exp 0)  ; maybe going denormalized       
288           (if (<= new-exp (- IEEE-double-float-digits))
289             0.0d0 ; should this be underflow? - should just be normal and result is fn of current fpu-mode
290             ;(error "Can't scale ~s by ~s." float int) ; should signal something                     
291             (let ((result (%make-dfloat)))
292               (%copy-double-float float result)
293               (set-%double-float-exp result 1) ; scale by float-exp -1
294               (%%scale-dfloat! result (+ IEEE-double-float-bias (+ float-exp int)) result)             
295               result))
296           (if (> new-exp IEEE-double-float-normal-exponent-max) 
297             (error (make-condition 'floating-point-overflow
298                                    :operation 'scale-float
299                                    :operands (list float int)))
300             (let ((new-float (%make-dfloat)))
301               (%copy-double-float float new-float)
302               (set-%double-float-exp new-float new-exp)
303               new-float))))))
304    (short-float
305     (let* ((float-exp (%short-float-exp float))
306            (new-exp (+ float-exp int)))
307       (if (eq 0 float-exp) ; already denormalized?
308         (if (%short-float-zerop float)
309           float
310           #+32-bit-target
311           (let ((result (%make-sfloat)))
312             (%%scale-sfloat! float (+ (1+ IEEE-single-float-bias) int) result))
313           #+64-bit-target
314           (%%scale-sfloat float (+ (1+ IEEE-single-float-bias) int)))
315         (if (<= new-exp 0)  ; maybe going denormalized       
316           (if (<= new-exp (- IEEE-single-float-digits))
317             ;; should this be underflow? - should just be normal and
318             ;; result is fn of current fpu-mode (error "Can't scale
319             ;; ~s by ~s." float int) ; should signal something
320             0.0s0
321             #+32-bit-target
322             (let ((result (%make-sfloat)))
323               (%copy-short-float float result)
324               (set-%short-float-exp result 1) ; scale by float-exp -1
325               (%%scale-sfloat! result (+ IEEE-single-float-bias (+ float-exp int)) result)             
326               result)
327             #+64-bit-target
328             (%%scale-sfloat (set-%short-float-exp float 1)
329                             (+ IEEE-single-float-bias (+ float-exp int))))
330           (if (> new-exp IEEE-single-float-normal-exponent-max) 
331             (error (make-condition 'floating-point-overflow
332                                    :operation 'scale-float
333                                    :operands (list float int)))
334             #+32-bit-target
335             (let ((new-float (%make-sfloat)))
336               (%copy-short-float float new-float)
337               (set-%short-float-exp new-float new-exp)
338               new-float)
339             #+64-bit-target
340             (set-%short-float-exp float new-exp))))))))
341
342(defun %copy-float (f)
343  ;Returns a freshly consed float.  float can also be a macptr.
344  (cond ((double-float-p f) (%copy-double-float f (%make-dfloat)))
345        ((macptrp f)
346         (let ((float (%make-dfloat)))
347           (%copy-ptr-to-ivector f 0 float (* 4 target::double-float.value-cell) 8)
348           float))
349        (t (error "Illegal arg ~s to %copy-float" f))))
350
351(defun float-precision (float)     ; not used - not in cltl2 index ?
352  "Return a non-negative number of significant digits in its float argument.
353  Will be less than FLOAT-DIGITS if denormalized or zero."
354  (number-case float
355     (double-float
356      (if (eq 0 (%double-float-exp float))
357        (if (not (%double-float-zerop float))
358        ; denormalized
359          (- IEEE-double-float-mantissa-width (dfloat-significand-zeros float))
360          0)
361        IEEE-double-float-digits))
362     (short-float 
363      (if (eq 0 (%short-float-exp float))
364        (if (not (%short-float-zerop float))
365        ; denormalized
366          (- IEEE-single-float-mantissa-width (sfloat-significand-zeros float))
367          0)
368        IEEE-single-float-digits))))
369
370
371(defun %double-float (number &optional result)
372  ;(require-null-or-double-float-sym result)
373  ; use number-case when macro is common
374  (number-case number
375    (double-float
376     (if result 
377       (%copy-double-float number result)
378         number))
379    (short-float
380     (%short-float->double-float number (or result (%make-dfloat))))
381    (fixnum
382     (%fixnum-dfloat number (or result (%make-dfloat))))
383    (bignum (%bignum-dfloat number result))
384    (ratio 
385     (if (not result)(setq result (%make-dfloat)))
386     (let* ((num (%numerator number))
387            (den (%denominator number)))
388       ; dont error if result is floatable when either top or bottom is not.
389       ; maybe do usual first, catching error
390       (if (not (or (bignump num)(bignump den)))
391         (with-stack-double-floats ((fnum num)
392                                        (fden den))       
393             (%double-float/-2! fnum fden result))
394         (let* ((numlen (integer-length num))
395                (denlen (integer-length den))
396                (exp (- numlen denlen))
397                (minusp (minusp num)))
398           (if (and (<= numlen IEEE-double-float-bias)
399                    (<= denlen IEEE-double-float-bias)
400                    #|(not (minusp exp))|# 
401                    (<= (abs exp) IEEE-double-float-mantissa-width))
402             (with-stack-double-floats ((fnum num)
403                                            (fden den))
404       
405               (%double-float/-2! fnum fden result))
406             (if (> exp IEEE-double-float-mantissa-width)
407               (progn  (%double-float (round num den) result))               
408               (if (>= exp 0)
409                 ; exp between 0 and 53 and nums big
410                 (let* ((shift (- IEEE-double-float-digits exp))
411                        (num (if minusp (- num) num))
412                        (int (round (ash num shift) den)) ; gaak
413                        (intlen (integer-length int))
414                        (new-exp (+ intlen (- IEEE-double-float-bias shift))))
415                   
416                   (when (> intlen IEEE-double-float-digits)
417                     (setq shift (1- shift))
418                     (setq int (round (ash num shift) den))
419                     (setq intlen (integer-length int))
420                     (setq new-exp (+ intlen (- IEEE-double-float-bias shift))))
421                   (when (> new-exp 2046)
422                     (error (make-condition 'floating-point-overflow
423                                            :operation 'double-float
424                                            :operands (list number))))
425                   (make-float-from-fixnums (ldb (byte 25 (- intlen 25)) int)
426                                            (ldb (byte 28 (max (- intlen 53) 0)) int)
427                                            new-exp ;(+ intlen (- IEEE-double-float-bias 53))
428                                            (if minusp -1 1)
429                                            result))
430                 ; den > num - exp negative
431                 (progn 
432                   (float-rat-neg-exp num den (if minusp -1 1) result)))))))))))
433
434
435#+32-bit-target
436(defun %short-float-ratio (number &optional result)
437  (if (not result)(setq result (%make-sfloat)))
438  (let* ((num (%numerator number))
439         (den (%denominator number)))
440    ;; dont error if result is floatable when either top or bottom is
441    ;; not.  maybe do usual first, catching error
442    (if (not (or (bignump num)(bignump den)))
443      (target::with-stack-short-floats ((fnum num)
444                                       (fden den))       
445        (%short-float/-2! fnum fden result))
446      (let* ((numlen (integer-length num))
447             (denlen (integer-length den))
448             (exp (- numlen denlen))
449             (minusp (minusp num)))
450        (if (and (<= numlen IEEE-single-float-bias)
451                 (<= denlen IEEE-single-float-bias)
452                 #|(not (minusp exp))|# 
453                 (<= (abs exp) IEEE-single-float-mantissa-width))
454          (target::with-stack-short-floats ((fnum num)
455                                           (fden den))
456            (%short-float/-2! fnum fden result))
457          (if (> exp IEEE-single-float-mantissa-width)
458            (progn  (%short-float (round num den) result))               
459            (if (>= exp 0)
460              ; exp between 0 and 23 and nums big
461              (let* ((shift (- IEEE-single-float-digits exp))
462                     (num (if minusp (- num) num))
463                     (int (round (ash num shift) den)) ; gaak
464                     (intlen (integer-length int))
465                     (new-exp (+ intlen (- IEEE-single-float-bias shift))))
466                (when (> intlen IEEE-single-float-digits)
467                  (setq shift (1- shift))
468                  (setq int (round (ash num shift) den))
469                  (setq intlen (integer-length int))
470                  (setq new-exp (+ intlen (- IEEE-single-float-bias shift))))
471                (when (> new-exp IEEE-single-float-normal-exponent-max)
472                  (error (make-condition 'floating-point-overflow
473                                         :operation 'short-float
474                                         :operands (list number))))
475                (make-short-float-from-fixnums 
476                   (ldb (byte IEEE-single-float-digits  (- intlen  IEEE-single-float-digits)) int)
477                   new-exp
478                   (if minusp -1 1)
479                   result))
480              ; den > num - exp negative
481              (progn 
482                (float-rat-neg-exp num den (if minusp -1 1) result t)))))))))
483
484#+64-bit-target
485(defun %short-float-ratio (number)
486  (let* ((num (%numerator number))
487         (den (%denominator number)))
488    ;; dont error if result is floatable when either top or bottom is
489    ;; not.  maybe do usual first, catching error
490    (if (not (or (bignump num)(bignump den)))
491      (/ (the short-float (%short-float num))
492         (the short-float (%short-float den)))
493      (let* ((numlen (integer-length num))
494             (denlen (integer-length den))
495             (exp (- numlen denlen))
496             (minusp (minusp num)))
497        (if (and (<= numlen IEEE-single-float-bias)
498                 (<= denlen IEEE-single-float-bias)
499                 #|(not (minusp exp))|# 
500                 (<= (abs exp) IEEE-single-float-mantissa-width))
501          (/ (the short-float (%short-float num))
502             (the short-float (%short-float den)))
503          (if (> exp IEEE-single-float-mantissa-width)
504            (progn  (%short-float (round num den)))
505            (if (>= exp 0)
506              ; exp between 0 and 23 and nums big
507              (let* ((shift (- IEEE-single-float-digits exp))
508                     (num (if minusp (- num) num))
509                     (int (round (ash num shift) den)) ; gaak
510                     (intlen (integer-length int))
511                     (new-exp (+ intlen (- IEEE-single-float-bias shift))))
512                (when (> intlen IEEE-single-float-digits)
513                  (setq shift (1- shift))
514                  (setq int (round (ash num shift) den))
515                  (setq intlen (integer-length int))
516                  (setq new-exp (+ intlen (- IEEE-single-float-bias shift))))
517                (when (> new-exp IEEE-single-float-normal-exponent-max)
518                  (error (make-condition 'floating-point-overflow
519                                         :operation 'short-float
520                                         :operands (list number))))
521                (make-short-float-from-fixnums 
522                   (ldb (byte IEEE-single-float-digits  (- intlen  IEEE-single-float-digits)) int)
523                   new-exp
524                   (if minusp 1 0)))
525              ; den > num - exp negative
526              (progn 
527                (float-rat-neg-exp num den (if minusp -1 1) nil t)))))))))
528
529
530#+32-bit-target
531(defun %short-float (number &optional result)
532  (number-case number
533    (short-float
534     (if result (%copy-short-float number result) number))
535    (double-float
536     (%double-float->short-float number (or result (%make-sfloat))))
537    (fixnum
538     (%fixnum-sfloat number (or result (%make-sfloat))))
539    (bignum
540     (%bignum-sfloat number (or result (%make-sfloat))))
541    (ratio
542     (%short-float-ratio number result))))
543
544#+64-bit-target
545(defun %short-float (number)
546  (number-case number
547    (short-float number)
548    (double-float (%double-float->short-float number))
549    (fixnum (%fixnum-sfloat number))
550    (bignum (%bignum-sfloat number))
551    (ratio (%short-float-ratio number))))
552
553
554(defun float-rat-neg-exp (integer divisor sign &optional result short)
555  (if (minusp sign)(setq integer (- integer)))       
556  (let* ((integer-length (integer-length integer))
557         ;; make sure we will have enough bits in the quotient
558         ;; (and a couple extra for rounding)
559         (shift-factor (+ (- (integer-length divisor) integer-length) (if short 28 60))) ; fix
560         (scaled-integer integer))
561    (if (plusp shift-factor)
562      (setq scaled-integer (ash integer shift-factor))
563      (setq divisor (ash divisor (- shift-factor)))  ; assume div > num
564      )
565    ;(pprint (list shift-factor scaled-integer divisor))
566    (multiple-value-bind (quotient remainder)(floor scaled-integer divisor)
567      (unless (zerop remainder) ; whats this - tells us there's junk below
568        (setq quotient (logior quotient 1)))
569      ;; why do it return 2 values?
570      (values (float-and-scale-and-round sign quotient (- shift-factor)  short result)))))
571
572
573
574;;; when is (negate-bignum (bignum-ashift-right big)) ; can't negate
575;;; in place cause may get bigger cheaper than (negate-bignum big) - 6
576;;; 0r 8 digits ; 8 longs so win if digits > 7 or negate it on the
577;;; stack
578
579(defun %bignum-dfloat (big &optional result) 
580  (let* ((minusp (bignum-minusp big)))
581    (flet 
582      ((doit (new-big)
583         (let* ((int-len (bignum-integer-length new-big)))
584           (when (>= int-len (- 2047 IEEE-double-float-bias)) ; args?
585             (error (make-condition 'floating-point-overflow 
586                                    :operation 'float :operands (list big))))
587           (if (> int-len 53)
588             (let* ((hi (ldb (byte 25  (- int-len  25)) new-big))
589                    (lo (ldb (byte 28 (- int-len 53)) new-big)))
590               ;(print (list new-big hi lo))
591               (when (and (logbitp (- int-len 54) new-big)  ; round bit
592                          (or (%ilogbitp 0 lo)    ; oddp
593                              ;; or more bits below round
594                              (%i< (one-bignum-factor-of-two new-big) (- int-len 54))))
595                 (if (eq lo #xfffffff)
596                   (setq hi (1+ hi) lo 0)
597                   (setq lo (1+ lo)))
598                 (when (%ilogbitp 25 hi) ; got bigger
599                   (setq int-len (1+ int-len))
600                   (let ((bit (%ilogbitp 0 hi)))
601                     (setq hi (%ilsr 1 hi))
602                     (setq lo (%ilsr 1 lo))
603                     (if bit (setq lo (%ilogior #x8000000 lo))))))
604               (make-float-from-fixnums hi lo (+ IEEE-double-float-bias int-len)(if minusp -1 1) result))
605             (let* ((hi (ldb (byte 25  (- int-len  25)) new-big))
606                    (lobits (min (- int-len 25) 28))
607                    (lo (ldb (byte lobits (- int-len (+ lobits 25))) new-big)))
608               (if (< lobits 28) (setq lo (ash lo (- 28 lobits))))
609               (make-float-from-fixnums hi lo (+ IEEE-double-float-bias int-len) (if minusp -1 1) result))))))
610      (declare (dynamic-extent #'doit))
611      (with-one-negated-bignum-buffer big doit))))
612
613#+32-bit-target
614(defun %bignum-sfloat (big &optional result) 
615  (let* ((minusp (bignum-minusp big)))
616    (flet 
617      ((doit (new-big)
618         (let* ((int-len (bignum-integer-length new-big)))
619           (when (>= int-len (- 255 IEEE-single-float-bias)) ; args?
620             (error (make-condition 'floating-point-overflow 
621                                    :operation 'float :operands (list big 1.0s0))))
622           (if t ;(> int-len IEEE-single-float-digits) ; always true
623             (let* ((lo (ldb (byte IEEE-single-float-digits  (- int-len  IEEE-single-float-digits)) new-big)))
624               (when (and (logbitp (- int-len 25) new-big)  ; round bit
625                          (or (%ilogbitp 0 lo)    ; oddp
626                              ; or more bits below round
627                              (%i< (one-bignum-factor-of-two new-big) (- int-len 25))))
628                 (setq lo (1+ lo))
629                 (when (%ilogbitp 24 lo) ; got bigger
630                   (setq int-len (1+ int-len))
631                   (setq lo (%ilsr 1 lo))))
632               (make-short-float-from-fixnums  lo (+ IEEE-single-float-bias int-len)(if minusp -1 1) result))
633             ))))
634      (declare (dynamic-extent #'doit))
635      (with-one-negated-bignum-buffer big doit))))
636
637
638#+64-bit-target
639(defun %bignum-sfloat (big) 
640  (let* ((minusp (bignum-minusp big)))
641    (flet 
642      ((doit (new-big)
643         (let* ((int-len (bignum-integer-length new-big)))
644           (when (>= int-len (- 255 IEEE-single-float-bias)) ; args?
645             (error (make-condition 'floating-point-overflow 
646                                    :operation 'float :operands (list big 1.0s0))))
647           (if t ;(> int-len IEEE-single-float-digits) ; always true
648             (let* ((lo (ldb (byte IEEE-single-float-digits  (- int-len  IEEE-single-float-digits)) new-big)))
649               (when (and (logbitp (- int-len 25) new-big)  ; round bit
650                          (or (%ilogbitp 0 lo)    ; oddp
651                              ; or more bits below round
652                              (%i< (one-bignum-factor-of-two new-big) (- int-len 25))))
653                 (setq lo (1+ lo))
654                 (when (%ilogbitp 24 lo) ; got bigger
655                   (setq int-len (1+ int-len))
656                   (setq lo (%ilsr 1 lo))))
657               (make-short-float-from-fixnums  lo (+ IEEE-single-float-bias int-len)(if minusp -1 1)))
658             ))))
659      (declare (dynamic-extent #'doit))
660      (with-one-negated-bignum-buffer big doit))))
661
662
663
664
665(defun %fixnum-dfloat (fix &optional result) 
666  (if (eq 0 fix) 
667    (if result (%copy-double-float 0.0d0 result) 0.0d0)
668    (progn
669      (when (not result)(setq result (%make-dfloat)))
670      ; it better return result
671      (%int-to-dfloat fix result))))
672
673
674#+32-bit-target
675(defun %fixnum-sfloat (fix &optional result)
676  (if (eq 0 fix)
677    (if result (%copy-short-float 0.0s0 result) 0.0s0)
678    (%int-to-sfloat! fix (or result (%make-sfloat)))))
679
680#+64-bit-target
681(defun %fixnum-sfloat (fix)
682  (if (eq 0 fix)
683    0.0s0
684    (%int-to-sfloat fix)))
685
686;;; Transcendental functions.
687(defun sin (x)
688  "Return the sine of NUMBER."
689  (if (complexp x)
690    (let* ((r (realpart x))
691           (i (imagpart x)))
692      (complex (* (sin r) (cosh i))
693               (* (cos r) (sinh i))))
694    (if (typep x 'double-float)
695      (%double-float-sin! x (%make-dfloat))
696      #+32-bit-target
697      (target::with-stack-short-floats ((sx x))
698        (%single-float-sin! sx (%make-sfloat)))
699      #+64-bit-target
700      (%single-float-sin (%short-float x)))))
701
702(defun cos (x)
703  "Return the cosine of NUMBER."
704  (if (complexp x)
705    (let* ((r (realpart x))
706           (i (imagpart x)))
707      (complex (* (cos r) (cosh i))
708               (- (* (sin r) (sinh i)))))
709    (if (typep x 'double-float)
710      (%double-float-cos! x (%make-dfloat))
711      #+32-bit-target
712      (target::with-stack-short-floats ((sx x))
713        (%single-float-cos! sx (%make-sfloat)))
714      #+64-bit-target
715      (%single-float-cos (%short-float x)))))
716
717(defun tan (x)
718  "Return the tangent of NUMBER."
719  (if (complexp x)
720    (/ (sin x) (cos x))
721    (if (typep x 'double-float)
722      (%double-float-tan! x (%make-dfloat))
723      #+32-bit-target
724      (target::with-stack-short-floats ((sx x))
725        (%single-float-tan! sx (%make-sfloat)))
726      #+64-bit-target
727      (%single-float-tan (%short-float x))
728      )))
729
730
731
732
733(defun atan (y &optional (x nil x-p))
734  "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
735  (if x-p
736    (if (or (typep x 'double-float)
737            (typep y 'double-float))
738      (with-stack-double-floats ((dy y)
739                                 (dx x))
740        (%df-atan2 dy dx))
741      #+32-bit-target
742      (target::with-stack-short-floats ((sy y)
743                                (sx x))
744        (%sf-atan2! sy sx))
745      #+64-bit-target
746      (%sf-atan2 (%short-float y) (%short-float x)))
747    (if (typep y 'complex)
748      (let* ((iy (* (sqrt -1) y)))
749             (/ (- (log (+ 1 iy)) (log (- 1 iy)))
750                #c(0 2)))
751      (if (typep y 'double-float)
752        (%double-float-atan! y (%make-dfloat))
753        #+32-bit-target
754        (target::with-stack-short-floats ((sy y))
755          (%single-float-atan! sy (%make-sfloat)))
756        #+64-bit-target
757        (%single-float-atan (%short-float y))
758        ))))
759
760
761
762(defun log (x &optional (b nil b-p))
763  "Return the logarithm of NUMBER in the base BASE, which defaults to e."
764  (if b-p
765    (if (zerop b)
766      (if (zerop x)
767        (report-bad-arg x '(not (satisfies zerop) ))
768        (if (floatp x) (float 0.0d0 x) 0))
769      (/ (log-e x) (log-e b)))
770    (log-e x)))
771
772(defun log-e (x)
773  (cond 
774    ((bignump x)
775     (if (minusp x)
776       (complex (log-e (- x)) pi)
777       (let* ((base1 3)
778              (guess (floor (1- (integer-length x))
779                            (log base1 2)))
780              (guess1 (* guess (log-e base1))))
781         (+ guess1 (log-e (/ x (expt base1 guess)))))))
782    ((and (ratiop x) 
783          (or (> x most-positive-short-float)
784              (< x most-negative-short-float)))
785     (- (log-e (%numerator x)) (log-e (%denominator x))))
786    ((typep x 'complex)
787     (complex (log-e (abs x)) (phase x)))
788    ((typep x 'double-float)
789     (with-stack-double-floats ((dx x))
790       (if (minusp x)
791         (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi)
792         (%double-float-log! dx (%make-dfloat)))))
793    (t
794     #+32-bit-target
795     (target::with-stack-short-floats ((sx x))
796       (if (minusp x)
797         (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat))
798                  #.(coerce pi 'short-float))
799         (%single-float-log! sx (%make-sfloat))))
800     #+64-bit-target
801     (if (minusp x)
802       (complex (%single-float-log (%short-float-abs (%short-float x))) #.(coerce pi 'single-float))
803       (%single-float-log (%short-float x)))
804     )))
805
806
807
808(defun exp (x)
809  "Return e raised to the power NUMBER."
810  (typecase x
811    (complex (* (exp (realpart x)) (cis (imagpart x))))
812    (double-float (%double-float-exp! x (%make-dfloat)))
813    (t
814     #+32-bit-target
815     (target::with-stack-short-floats ((sx x))
816       (%single-float-exp! sx (%make-sfloat)))
817     #+64-bit-target
818     (%single-float-exp (%short-float x)))))
819
820
821(defun positive-realpart-p (n)
822  (> (realpart n) 0))
823
824(defun expt (b e)
825  "Return BASE raised to the POWER."
826  (cond ((zerop e) (1+ (* b e)))
827        ((integerp e)
828         (if (minusp e) (/ 1 (%integer-power b (- e))) (%integer-power b e)))
829        ((zerop b)
830         (if (plusp (realpart e)) b (report-bad-arg e '(satisfies positive-realpart-p))))
831        ((and (realp b) (plusp b) (realp e))
832         (if (or (typep b 'double-float)
833                 (typep e 'double-float))
834           (with-stack-double-floats ((b1 b)
835                                      (e1 e))
836             (%double-float-expt! b1 e1 (%make-dfloat)))
837           #+32-bit-target
838           (target::with-stack-short-floats ((b1 b)
839                                     (e1 e))
840             (%single-float-expt! b1 e1 (%make-sfloat)))
841           #+64-bit-target
842           (%single-float-expt (%short-float b) (%short-float e))
843           ))
844        ((typep (realpart e) 'double-float)
845         ;; Avoid intermediate single-float result from LOG
846         (let ((promoted-base (* 1d0 b)))
847           (exp (* e (log promoted-base)))))
848        (t (exp (* e (log b))))))
849
850
851
852(defun sqrt (x &aux a b)
853  "Return the square root of NUMBER."
854  (cond ((zerop x) x)
855        ((complexp x) (* (sqrt (abs x)) (cis (/ (phase x) 2))))         
856        ((minusp x) (complex 0 (sqrt (- x))))
857        ((floatp x)
858         (fsqrt x))
859        ((and (integerp x) (eql x (* (setq a (isqrt x)) a))) a)
860        ((and (ratiop x)
861              (let ((n (numerator x))
862                    d)
863                (and (eql n (* (setq a (isqrt n)) a))
864                     (eql (setq d (denominator x))
865                          (* (setq b (isqrt d)) b)))))
866         (/ a b))         
867        (t
868         #+32-bit-target
869         (target::with-stack-short-floats ((f1))
870           (fsqrt (%short-float x f1)))
871         #+64-bit-target
872         (fsqrt (%short-float x)))))
873
874
875
876(defun asin (x)
877  "Return the arc sine of NUMBER."
878  (number-case x
879    (complex
880      (let ((sqrt-1-x (sqrt (- 1 x)))
881            (sqrt-1+x (sqrt (+ 1 x))))
882        (complex (atan (/ (realpart x)
883                          (realpart (* sqrt-1-x sqrt-1+x))))
884                 (asinh (imagpart (* (conjugate sqrt-1-x)
885                                     sqrt-1+x))))))
886    (double-float
887     (locally (declare (type double-float x))
888       (if (and (<= -1.0d0 x)
889                (<= x 1.0d0))
890         (%double-float-asin! x (%make-dfloat))
891         (let* ((temp (+ (complex -0.0d0 x)
892                         (sqrt (- 1.0d0 (the double-float (* x x)))))))
893           (complex (phase temp) (- (log (abs temp))))))))
894    ((short-float rational)
895     #+32-bit-target
896     (let* ((x1 (%make-sfloat)))
897       (declare (dynamic-extent x1))
898       (if (and (realp x) 
899                (<= -1.0s0 (setq x (%short-float x x1)))
900                (<= x 1.0s0))
901         (%single-float-asin! x1 (%make-sfloat))
902         (progn
903           (setq x (+ (complex (- (imagpart x)) (realpart x))
904                      (sqrt (- 1 (* x x)))))
905           (complex (phase x) (- (log (abs x)))))))
906     #+64-bit-target
907     (if (and (realp x) 
908              (<= -1.0s0 (setq x (%short-float x)))
909              (<= x 1.0s0))
910         (%single-float-asin x)
911         (progn
912           (setq x (+ (complex (- (imagpart x)) (realpart x))
913                      (sqrt (- 1 (* x x)))))
914           (complex (phase x) (- (log (abs x))))))
915     )))
916
917
918(eval-when (:execute :compile-toplevel)
919  (defconstant single-float-pi (coerce pi 'single-float))
920  (defconstant double-float-half-pi (asin 1.0d0))
921  (defconstant single-float-half-pi (asin 1.0f0))
922)
923
924
925
926(defun acos (x)
927  "Return the arc cosine of NUMBER."
928  (number-case x
929    (complex
930     (let ((sqrt-1+x (sqrt (+ 1 x)))
931           (sqrt-1-x (sqrt (- 1 x))))
932       (complex (* 2 (atan (/ (realpart sqrt-1-x)
933                              (realpart sqrt-1+x))))
934                (asinh (imagpart (* (conjugate sqrt-1+x)
935                                    sqrt-1-x))))))
936   
937    (double-float
938     (locally (declare (type double-float x))
939       (if (and (<= -1.0d0 x)
940                (<= x 1.0d0))
941         (%double-float-acos! x (%make-dfloat))
942         (- double-float-half-pi (asin x)))))
943    ((short-float rational)
944     #+32-bit-target
945     (target::with-stack-short-floats ((sx x))
946        (locally
947            (declare (type short-float sx))
948          (if (and (<= -1.0s0 sx)
949                   (<= sx 1.0s0))
950            (%single-float-acos! sx (%make-sfloat))
951            (- single-float-half-pi (asin sx)))))
952     #+64-bit-target
953     (let* ((sx (%short-float x)))
954       (declare (type short-float sx))
955       (if (and (<= -1.0s0 sx)
956                (<= sx 1.0s0))
957         (%single-float-acos sx)
958         (- single-float-half-pi (asin sx))))
959     )))
960
961
962(defun fsqrt (x)
963  (etypecase x
964    (double-float (%double-float-sqrt! x (%make-dfloat)))
965    (single-float
966     #+32-bit-target
967     (%single-float-sqrt! x (%make-sfloat))
968     #+64-bit-target
969     (%single-float-sqrt x))))
970
971
972
973(defun %df-atan2 (y x)
974  (if (zerop x)
975    (if (zerop y)
976      (if (plusp (float-sign x))
977        (if (eql y -0.0d0)
978          -0.0d0
979          0.0d0)
980        (float-sign y pi))
981      (float-sign y double-float-half-pi))
982    (%double-float-atan2! y x (%make-dfloat))))
983
984#+32-bit-target
985(defun %sf-atan2! (y x)
986  (if (zerop x)
987    (if (zerop y)
988      (if (plusp (float-sign x))
989        ;; Don't return Y (which may be stack-consed) here.
990        ;; We know that (ZEROP Y) is true, so:
991        (if (eql y -0.0s0)
992          -0.0s0
993          0.0s0)
994        (float-sign y single-float-pi))
995      (float-sign y single-float-half-pi))
996    (%single-float-atan2! y x (%make-sfloat))))
997
998#+64-bit-target
999(defun %sf-atan2 (y x)
1000  (if (zerop x)
1001    (if (zerop y)
1002      (if (plusp (float-sign x))
1003        y
1004        (float-sign y single-float-pi))
1005      (float-sign y single-float-half-pi))
1006    (%single-float-atan2 y x)))
1007
1008#+64-bit-target
1009(defun %short-float-exp (n)
1010  (let* ((bits (single-float-bits n)))
1011    (declare (type (unsigned-byte 32) bits))
1012    (ldb (byte IEEE-single-float-exponent-width IEEE-single-float-exponent-offset) bits)))
1013
1014
1015#+64-bit-target
1016(defun set-%short-float-exp (float exp)
1017  (host-single-float-from-unsigned-byte-32
1018   (dpb exp
1019        (byte IEEE-single-float-exponent-width
1020              IEEE-single-float-exponent-offset)
1021        (the (unsigned-byte 32) (single-float-bits float)))))
1022
1023#+64-bit-target
1024(defun %%scale-sfloat (float int)
1025  (* (the single-float float)
1026     (the single-float (host-single-float-from-unsigned-byte-32
1027                        (dpb int
1028                             (byte IEEE-single-float-exponent-width
1029                                   IEEE-single-float-exponent-offset)
1030                             0)))))
1031
1032#+64-bit-target
1033(defun %double-float-exp (n)
1034  (let* ((highword (double-float-bits n)))
1035    (declare (fixnum highword))
1036    (logand (1- (ash 1 IEEE-double-float-exponent-width))
1037            (ash highword (- (- IEEE-double-float-exponent-offset 32))))))
1038
1039#+64-bit-target
1040(defun set-%double-float-exp (float exp)
1041  (let* ((highword (double-float-bits float)))
1042    (declare (fixnum highword))
1043    (setf (uvref float target::double-float.val-high-cell)
1044          (dpb exp
1045               (byte IEEE-double-float-exponent-width
1046                     (- IEEE-double-float-exponent-offset 32))
1047               highword))
1048    exp))
1049
1050#+64-bit-target
1051(defun %integer-decode-double-float (f)
1052  (multiple-value-bind (hiword loword) (double-float-bits f)
1053    (declare (type (unsigned-byte 32) hiword loword))
1054    (let* ((exp (ldb (byte IEEE-double-float-exponent-width
1055                           (- IEEE-double-float-exponent-offset 32))
1056                     hiword))
1057           (mantissa (logior
1058                      (the fixnum
1059                        (dpb (ldb (byte (- IEEE-double-float-mantissa-width 32)
1060                                        IEEE-double-float-mantissa-offset)
1061                                  hiword)
1062                             (byte (- IEEE-double-float-mantissa-width 32)
1063                                   32)
1064                             loword))
1065                      (if (zerop exp)
1066                        0
1067                        (ash 1 IEEE-double-float-hidden-bit))))
1068           (sign (if (logbitp 31 hiword) -1 1)))
1069      (declare (fixnum exp mantissa sign))
1070      (values (ldb (byte 25 28) mantissa)
1071              (ldb (byte 28 0) mantissa)
1072              exp
1073              sign))))
1074
1075;;; end of l0-float.lisp
Note: See TracBrowser for help on using the repository browser.