source: release/1.2/source/level-0/l0-float.lisp @ 10208

Last change on this file since 10208 was 10208, checked in by rme, 11 years ago

Merge r10013 (bug in complex cos).

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