source: release/1.9/source/level-0/l0-float.lisp @ 15640

Last change on this file since 15640 was 15640, checked in by gb, 6 years ago

Avoid duplicate/different definitions of some constants.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.7 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  (defconstant single-float-pi (coerce pi 'single-float))
26  (defconstant double-float-half-pi (/ pi 2))
27  (defconstant single-float-half-pi (coerce (/ pi 2) 'single-float))
28  (defconstant single-float-log2 0.6931472)                ; (log 2)
29  (defconstant double-float-log2 0.6931471805599453d0)     ; (log 2.0d0)
30  (defconstant double-float-log2^23 15.942385152878742d0)  ; (log (expt 2 23))
31)
32
33;;; used by float reader
34(defun make-float-from-fixnums (hi lo exp sign &optional result)
35  ;;(require-null-or-double-float-sym result)
36  ;; maybe nuke all these require-types?
37  ;;(setq hi (require-type hi 'fixnum))
38  ;;(setq lo (require-type lo 'fixnum))
39  ;;(setq exp (require-type exp 'fixnum))
40  ;;(setq sign (require-type sign 'fixnum))
41  (let ((the-float (or result (%make-dfloat))))
42    (%make-float-from-fixnums the-float hi lo exp sign)
43    the-float))
44
45
46#+32-bit-target
47(defun make-short-float-from-fixnums (significand biased-exp sign &optional result)
48  (%make-short-float-from-fixnums (or result (%make-sfloat)) significand biased-exp sign))
49
50#+64-bit-target
51(defun make-short-float-from-fixnums (significand biased-exp sign)
52  (declare (fixnum significand biased-exp sign))
53  (host-single-float-from-unsigned-byte-32
54   (logior
55    (the fixnum (if (< sign 0) (ash 1 31) 0))
56    (the fixnum (ash biased-exp IEEE-single-float-exponent-offset))
57    (the fixnum (logand significand
58                        (1- (ash 1 IEEE-single-float-hidden-bit)))))))
59
60
61(defun float-sign (n1 &optional n2) ; second arg silly
62  "Return a floating-point number that has the same sign as
63   FLOAT1 and, if FLOAT2 is given, has the same absolute value
64   as FLOAT2."
65  (if (and n2 (not (typep n2 'float)))
66    (setq n2 (require-type n2 'float)))
67  (number-case n1
68    (double-float                       
69     (if (%double-float-sign n1) 
70       (if n2
71         (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) n2 (- n2))
72         -1.0d0)
73       (if n2
74         (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) (- n2) n2)
75         1.0d0)))
76    (short-float
77     (if (%short-float-sign n1)
78       (if n2
79         (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) n2 (- n2))
80         -1.0s0)
81       (if n2
82         (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) (- n2) n2)
83         1.0s0)))))
84
85
86
87(defun %double-float-minusp (n)
88  (and (%double-float-sign n)(not (%double-float-zerop n))))
89
90(defun %short-float-minusp (n)
91  (and (%short-float-sign n) (not (%short-float-zerop n))))
92
93(defun %double-float-abs (n)
94  (if (not (%double-float-sign n))
95    n 
96    (%%double-float-abs! n (%make-dfloat))))
97
98#+32-bit-target
99(defun %short-float-abs (n)
100  (if (not (%short-float-sign n))
101    n 
102    (%%short-float-abs! n (%make-sfloat))))
103
104(defun fixnum-decode-float (n)
105  (etypecase n
106    (double-float (%integer-decode-double-float n))))
107
108(defun nan-or-infinity-p (n)
109  (etypecase n
110    (double-float (eq 2047 (%double-float-exp n)))
111    (short-float (eq 255 (%short-float-exp n)))))
112
113; not sure this is right
114(defun infinity-p (n)
115  (etypecase n
116    (double-float (multiple-value-bind (hi lo exp)(fixnum-decode-float n)
117                    (and (eq 2047 exp)
118                         (eq #x1000000 hi)
119                         (eq 0 lo))))
120    (short-float
121     #+32-bit-target
122     (multiple-value-bind (high low)(%sfloat-hwords n)
123                  (let*  ((mantissa (%ilogior2 low (%ilsl 16 (%ilogand2 high #x007F))))
124                          (exp (%ilsr 7 (%ilogand2 high #x7F80))))
125                    (and (eq exp 255)
126                         (eq 0 mantissa))))
127     #+64-bit-target
128     (let* ((bits (single-float-bits n))
129            (exp (ldb (byte IEEE-single-float-exponent-width
130                            IEEE-single-float-exponent-offset)
131                      bits))
132            (mantissa (ldb (byte IEEE-single-float-mantissa-width
133                            IEEE-single-float-mantissa-offset)
134                           bits)))
135       (declare (fixnum bits exp mantissa))
136       (and (= exp 255)
137            (zerop mantissa))))))
138
139#+32-bit-target
140(defun fixnum-decode-short-float (float)
141  (multiple-value-bind (high low)(%sfloat-hwords float)
142    (let*  ((mantissa (%ilogior2 low (%ilsl 16 (%ilogand2 high #x007F))))
143            (exp (%ilsr 7 (%ilogand2 high #x7F80))))
144      (if (and (neq exp 0)(neq exp 255))(setq mantissa (%ilogior mantissa #x800000)))
145      (values mantissa exp (%ilsr 15 high)))))
146
147#+64-bit-target
148(defun fixnum-decode-short-float (float)
149  (let* ((bits (single-float-bits float)))
150    (declare (fixnum bits))
151    (let* ((mantissa (ldb (byte IEEE-single-float-mantissa-width
152                                IEEE-single-float-mantissa-offset)
153                          bits))
154           (exp (ldb (byte IEEE-single-float-exponent-width
155                           IEEE-single-float-exponent-offset)
156                     bits))
157           (sign (lsh bits -31)))
158      (declare (fixnum mantissa exp sign))
159      (unless (or (= exp 0) (= exp 255))
160        (setq mantissa (logior mantissa (ash 1 IEEE-single-float-hidden-bit))))
161      (values mantissa exp sign))))
162                 
163                   
164
165#+32-bit-target
166(defun integer-decode-double-float (n)
167  (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n)
168    ; is only 53 bits and positive so should be easy
169    ;(values (logior (ash hi 28) lo) exp sign)))
170    ; if denormalized, may fit in a fixnum
171    (setq exp (- exp (if (< hi #x1000000) 
172                       (+ IEEE-double-float-mantissa-width IEEE-double-float-bias)
173                       (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias)))))
174    (if (< hi (ash 1 (1- target::fixnumshift))) ; aka 2
175      (values (logior (ash hi 28) lo) exp sign)
176      ; might fit in 1 word?
177      (let ((big (%alloc-misc 2 target::subtag-bignum)))
178        (make-big-53 hi lo big)
179        (if (< hi #x1000000) (%normalize-bignum big))
180        (values big exp sign)))))
181
182#+64-bit-target
183(defun integer-decode-double-float (n)
184  (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n)
185    (setq exp (- exp (if (< hi #x1000000) 
186                       (+ IEEE-double-float-mantissa-width IEEE-double-float-bias)
187                       (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias)))))
188    (values (logior (ash hi 28) lo) exp sign)))
189   
190
191;;; actually only called when magnitude bigger than a fixnum
192#+32-bit-target
193(defun %truncate-double-float (n)
194  (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n)
195    (if (< exp (1+ IEEE-double-float-bias)) ; this is false in practice
196      0
197      (progn
198        (setq exp (- exp (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias))))
199        (if (eq sign 1)  ; positive
200          (logior (ash hi (+ 28 exp))(ash lo exp))
201          (if (<= exp 0) ; exp positive - negate before shift - else after
202            (let ((poo (logior (ash hi (+ 28 exp))(ash lo exp))))
203              (- poo))
204            (let ((poo (logior (ash hi 28) lo)))
205              (ash (- poo) exp))))))))
206
207#+64-bit-target
208(defun %truncate-double-float (n)
209  (multiple-value-bind (mantissa exp sign) (integer-decode-float n)
210    (* sign (ash mantissa exp))))
211
212
213
214; actually only called when bigger than a fixnum
215(defun %truncate-short-float (n)
216  (multiple-value-bind (mantissa exp sign)(fixnum-decode-short-float n)
217    (if (< exp (1+ IEEE-single-float-bias)) ; is magnitude less than 1 - false in practice
218      0
219      (progn
220        (setq exp (- exp (+ IEEE-single-float-mantissa-width (1+ IEEE-single-float-bias))))
221        (ash (if (eq sign 0) mantissa (- mantissa)) exp)))))
222
223(defun decode-float (n)
224  "Return three values:
225   1) a floating-point number representing the significand. This is always
226      between 0.5 (inclusive) and 1.0 (exclusive).
227   2) an integer representing the exponent.
228   3) -1.0 or 1.0 (i.e. the sign of the argument.)"
229  (number-case n
230    (double-float
231     (let* ((old-exp (%double-float-exp n))
232            (sign (if (%double-float-sign n) -1.0d0 1.0d0)))   
233       (if (eq 0 old-exp)
234         (if (%double-float-zerop n)
235           (values 0.0d0 0 sign)
236           (let* ((val (%make-dfloat))
237                  (zeros (dfloat-significand-zeros n)))
238             (%%double-float-abs! n val)
239             (%%scale-dfloat! val (+ 2 IEEE-double-float-bias zeros) val) ; get it normalized
240             (set-%double-float-exp val IEEE-double-float-bias)      ; then bash exponent
241             (values val (- old-exp zeros IEEE-double-float-bias) sign)))
242         (if (> old-exp IEEE-double-float-normal-exponent-max)
243           (error "Can't decode NAN or infinity ~s" n)
244           (let ((val (%make-dfloat)))
245             (%%double-float-abs! n val)
246             (set-%double-float-exp val IEEE-double-float-bias)
247             (values val (- old-exp IEEE-double-float-bias) sign))))))
248    (short-float
249     (let* ((old-exp (%short-float-exp n))
250            (sign (if (%short-float-sign n) -1.0s0 1.0s0)))
251       (if (eq 0 old-exp)
252         (if (%short-float-zerop n)
253           (values 0.0s0 0 sign)
254           #+32-bit-target
255           (let* ((val (%make-sfloat))
256                  (zeros (sfloat-significand-zeros n)))
257             (%%short-float-abs! n val)
258             (%%scale-sfloat! val (+ 2 IEEE-single-float-bias zeros) val) ; get it normalized
259             (set-%short-float-exp val IEEE-single-float-bias)      ; then bash exponent
260             (values val (- old-exp zeros IEEE-single-float-bias) sign))
261           #+64-bit-target
262           (let* ((zeros (sfloat-significand-zeros n))
263                  (val (%%scale-sfloat (%short-float-abs n)
264                                       (+ 2 IEEE-single-float-bias zeros))))
265             (values (set-%short-float-exp val IEEE-single-float-bias)
266                     (- old-exp zeros IEEE-single-float-bias) sign)))
267         (if (> old-exp IEEE-single-float-normal-exponent-max)
268           (error "Can't decode NAN or infinity ~s" n)
269           #+32-bit-target
270           (let ((val (%make-sfloat)))
271             (%%short-float-abs! n val)
272             (set-%short-float-exp val IEEE-single-float-bias)
273             (values val (- old-exp IEEE-single-float-bias) sign))
274           #+64-bit-target
275           (values (set-%short-float-exp (%short-float-abs n)
276                                         IEEE-single-float-bias)
277                   (- old-exp IEEE-single-float-bias) sign)))))))
278
279; (* float (expt 2 int))
280(defun scale-float (float int)
281  "Return the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
282  of precision or overflow."
283  (unless (fixnump int)(setq int (require-type int 'fixnum)))
284  (number-case float
285    (double-float
286     (let* ((float-exp (%double-float-exp float))
287            (new-exp (+ float-exp int)))
288       (if (eq 0 float-exp) ; already denormalized?
289         (if (%double-float-zerop float)
290           float 
291           (let ((result (%make-dfloat)))
292             (%%scale-dfloat! float (+ (1+ IEEE-double-float-bias) int) result)))
293         (if (<= new-exp 0)  ; maybe going denormalized       
294           (if (<= new-exp (- IEEE-double-float-digits))
295             0.0d0 ; should this be underflow? - should just be normal and result is fn of current fpu-mode
296             ;(error "Can't scale ~s by ~s." float int) ; should signal something                     
297             (let ((result (%make-dfloat)))
298               (%copy-double-float float result)
299               (set-%double-float-exp result 1) ; scale by float-exp -1
300               (%%scale-dfloat! result (+ IEEE-double-float-bias (+ float-exp int)) result)             
301               result))
302           (if (> new-exp IEEE-double-float-normal-exponent-max) 
303             (error (make-condition 'floating-point-overflow
304                                    :operation 'scale-float
305                                    :operands (list float int)))
306             (let ((new-float (%make-dfloat)))
307               (%copy-double-float float new-float)
308               (set-%double-float-exp new-float new-exp)
309               new-float))))))
310    (short-float
311     (let* ((float-exp (%short-float-exp float))
312            (new-exp (+ float-exp int)))
313       (if (eq 0 float-exp) ; already denormalized?
314         (if (%short-float-zerop float)
315           float
316           #+32-bit-target
317           (let ((result (%make-sfloat)))
318             (%%scale-sfloat! float (+ (1+ IEEE-single-float-bias) int) result))
319           #+64-bit-target
320           (%%scale-sfloat float (+ (1+ IEEE-single-float-bias) int)))
321         (if (<= new-exp 0)  ; maybe going denormalized       
322           (if (<= new-exp (- IEEE-single-float-digits))
323             ;; should this be underflow? - should just be normal and
324             ;; result is fn of current fpu-mode (error "Can't scale
325             ;; ~s by ~s." float int) ; should signal something
326             0.0s0
327             #+32-bit-target
328             (let ((result (%make-sfloat)))
329               (%copy-short-float float result)
330               (set-%short-float-exp result 1) ; scale by float-exp -1
331               (%%scale-sfloat! result (+ IEEE-single-float-bias (+ float-exp int)) result)             
332               result)
333             #+64-bit-target
334             (%%scale-sfloat (set-%short-float-exp float 1)
335                             (+ IEEE-single-float-bias (+ float-exp int))))
336           (if (> new-exp IEEE-single-float-normal-exponent-max) 
337             (error (make-condition 'floating-point-overflow
338                                    :operation 'scale-float
339                                    :operands (list float int)))
340             #+32-bit-target
341             (let ((new-float (%make-sfloat)))
342               (%copy-short-float float new-float)
343               (set-%short-float-exp new-float new-exp)
344               new-float)
345             #+64-bit-target
346             (set-%short-float-exp float new-exp))))))))
347
348(defun %copy-float (f)
349  ;Returns a freshly consed float.  float can also be a macptr.
350  (cond ((double-float-p f) (%copy-double-float f (%make-dfloat)))
351        ((macptrp f)
352         (let ((float (%make-dfloat)))
353           (%copy-ptr-to-ivector f 0 float (* 4 target::double-float.value-cell) 8)
354           float))
355        (t (error "Illegal arg ~s to %copy-float" f))))
356
357(defun float-precision (float)     ; not used - not in cltl2 index ?
358  "Return a non-negative number of significant digits in its float argument.
359  Will be less than FLOAT-DIGITS if denormalized or zero."
360  (number-case float
361     (double-float
362      (if (eq 0 (%double-float-exp float))
363        (if (not (%double-float-zerop float))
364        ; denormalized
365          (- IEEE-double-float-mantissa-width (dfloat-significand-zeros float))
366          0)
367        IEEE-double-float-digits))
368     (short-float 
369      (if (eq 0 (%short-float-exp float))
370        (if (not (%short-float-zerop float))
371        ; denormalized
372          (- IEEE-single-float-mantissa-width (sfloat-significand-zeros float))
373          0)
374        IEEE-single-float-digits))))
375
376
377(defun %double-float (number &optional result)
378  ;(require-null-or-double-float-sym result)
379  ; use number-case when macro is common
380  (number-case number
381    (double-float
382     (if result 
383       (%copy-double-float number result)
384         number))
385    (short-float
386     (%short-float->double-float number (or result (%make-dfloat))))
387    (fixnum
388     (%fixnum-dfloat number (or result (%make-dfloat))))
389    (bignum (%bignum-dfloat number result))
390    (ratio 
391     (if (not result)(setq result (%make-dfloat)))
392     (let* ((num (%numerator number))
393            (den (%denominator number)))
394       ; dont error if result is floatable when either top or bottom is not.
395       ; maybe do usual first, catching error
396       (if (not (or (bignump num)(bignump den)))
397         (with-stack-double-floats ((fnum num)
398                                        (fden den))       
399             (%double-float/-2! fnum fden result))
400         (let* ((numlen (integer-length num))
401                (denlen (integer-length den))
402                (exp (- numlen denlen))
403                (minusp (minusp num)))
404           (if (and (<= numlen IEEE-double-float-bias)
405                    (<= denlen IEEE-double-float-bias)
406                    #|(not (minusp exp))|# 
407                    (<= (abs exp) IEEE-double-float-mantissa-width))
408             (with-stack-double-floats ((fnum num)
409                                            (fden den))
410       
411               (%double-float/-2! fnum fden result))
412             (if (> exp IEEE-double-float-mantissa-width)
413               (progn  (%double-float (round num den) result))               
414               (if (>= exp 0)
415                 ; exp between 0 and 53 and nums big
416                 (let* ((shift (- IEEE-double-float-digits exp))
417                        (num (if minusp (- num) num))
418                        (int (round (ash num shift) den)) ; gaak
419                        (intlen (integer-length int))
420                        (new-exp (+ intlen (- IEEE-double-float-bias shift))))
421                   
422                   (when (> intlen IEEE-double-float-digits)
423                     (setq shift (1- shift))
424                     (setq int (round (ash num shift) den))
425                     (setq intlen (integer-length int))
426                     (setq new-exp (+ intlen (- IEEE-double-float-bias shift))))
427                   (when (> new-exp 2046)
428                     (error (make-condition 'floating-point-overflow
429                                            :operation 'double-float
430                                            :operands (list number))))
431                   (make-float-from-fixnums (ldb (byte 25 (- intlen 25)) int)
432                                            (ldb (byte 28 (max (- intlen 53) 0)) int)
433                                            new-exp ;(+ intlen (- IEEE-double-float-bias 53))
434                                            (if minusp -1 1)
435                                            result))
436                 ; den > num - exp negative
437                 (progn 
438                   (float-rat-neg-exp num den (if minusp -1 1) result)))))))))))
439
440
441#+32-bit-target
442(defun %short-float-ratio (number &optional result)
443  (if (not result)(setq result (%make-sfloat)))
444  (let* ((num (%numerator number))
445         (den (%denominator number)))
446    ;; dont error if result is floatable when either top or bottom is
447    ;; not.  maybe do usual first, catching error
448    (if (not (or (bignump num)(bignump den)))
449      (target::with-stack-short-floats ((fnum num)
450                                       (fden den))       
451        (%short-float/-2! fnum fden result))
452      (let* ((numlen (integer-length num))
453             (denlen (integer-length den))
454             (exp (- numlen denlen))
455             (minusp (minusp num)))
456        (if (and (<= numlen IEEE-single-float-bias)
457                 (<= denlen IEEE-single-float-bias)
458                 #|(not (minusp exp))|# 
459                 (<= (abs exp) IEEE-single-float-mantissa-width))
460          (target::with-stack-short-floats ((fnum num)
461                                           (fden den))
462            (%short-float/-2! fnum fden result))
463          (if (> exp IEEE-single-float-mantissa-width)
464            (progn  (%short-float (round num den) result))               
465            (if (>= exp 0)
466              ; exp between 0 and 23 and nums big
467              (let* ((shift (- IEEE-single-float-digits exp))
468                     (num (if minusp (- num) num))
469                     (int (round (ash num shift) den)) ; gaak
470                     (intlen (integer-length int))
471                     (new-exp (+ intlen (- IEEE-single-float-bias shift))))
472                (when (> intlen IEEE-single-float-digits)
473                  (setq shift (1- shift))
474                  (setq int (round (ash num shift) den))
475                  (setq intlen (integer-length int))
476                  (setq new-exp (+ intlen (- IEEE-single-float-bias shift))))
477                (when (> new-exp IEEE-single-float-normal-exponent-max)
478                  (error (make-condition 'floating-point-overflow
479                                         :operation 'short-float
480                                         :operands (list number))))
481                (make-short-float-from-fixnums 
482                   (ldb (byte IEEE-single-float-digits  (- intlen  IEEE-single-float-digits)) int)
483                   new-exp
484                   (if minusp -1 1)
485                   result))
486              ; den > num - exp negative
487              (progn 
488                (float-rat-neg-exp num den (if minusp -1 1) result t)))))))))
489
490#+64-bit-target
491(defun %short-float-ratio (number)
492  (let* ((num (%numerator number))
493         (den (%denominator number)))
494    ;; dont error if result is floatable when either top or bottom is
495    ;; not.  maybe do usual first, catching error
496    (if (not (or (bignump num)(bignump den)))
497      (/ (the short-float (%short-float num))
498         (the short-float (%short-float den)))
499      (let* ((numlen (integer-length num))
500             (denlen (integer-length den))
501             (exp (- numlen denlen))
502             (minusp (minusp num)))
503        (if (and (<= numlen IEEE-single-float-bias)
504                 (<= denlen IEEE-single-float-bias)
505                 #|(not (minusp exp))|# 
506                 (<= (abs exp) IEEE-single-float-mantissa-width))
507          (/ (the short-float (%short-float num))
508             (the short-float (%short-float den)))
509          (if (> exp IEEE-single-float-mantissa-width)
510            (progn  (%short-float (round num den)))
511            (if (>= exp 0)
512              ; exp between 0 and 23 and nums big
513              (let* ((shift (- IEEE-single-float-digits exp))
514                     (num (if minusp (- num) num))
515                     (int (round (ash num shift) den)) ; gaak
516                     (intlen (integer-length int))
517                     (new-exp (+ intlen (- IEEE-single-float-bias shift))))
518                (when (> intlen IEEE-single-float-digits)
519                  (setq shift (1- shift))
520                  (setq int (round (ash num shift) den))
521                  (setq intlen (integer-length int))
522                  (setq new-exp (+ intlen (- IEEE-single-float-bias shift))))
523                (when (> new-exp IEEE-single-float-normal-exponent-max)
524                  (error (make-condition 'floating-point-overflow
525                                         :operation 'short-float
526                                         :operands (list number))))
527                (make-short-float-from-fixnums 
528                   (ldb (byte IEEE-single-float-digits  (- intlen  IEEE-single-float-digits)) int)
529                   new-exp
530                   (if minusp 1 0)))
531              ; den > num - exp negative
532              (progn 
533                (float-rat-neg-exp num den (if minusp -1 1) nil t)))))))))
534
535
536#+32-bit-target
537(defun %short-float (number &optional result)
538  (number-case number
539    (short-float
540     (if result (%copy-short-float number result) number))
541    (double-float
542     (%double-float->short-float number (or result (%make-sfloat))))
543    (fixnum
544     (%fixnum-sfloat number (or result (%make-sfloat))))
545    (bignum
546     (%bignum-sfloat number (or result (%make-sfloat))))
547    (ratio
548     (%short-float-ratio number result))))
549
550#+64-bit-target
551(defun %short-float (number)
552  (number-case number
553    (short-float number)
554    (double-float (%double-float->short-float number))
555    (fixnum (%fixnum-sfloat number))
556    (bignum (%bignum-sfloat number))
557    (ratio (%short-float-ratio number))))
558
559
560(defun float-rat-neg-exp (integer divisor sign &optional result short)
561  (if (minusp sign)(setq integer (- integer)))       
562  (let* ((integer-length (integer-length integer))
563         ;; make sure we will have enough bits in the quotient
564         ;; (and a couple extra for rounding)
565         (shift-factor (+ (- (integer-length divisor) integer-length) (if short 28 60))) ; fix
566         (scaled-integer integer))
567    (if (plusp shift-factor)
568      (setq scaled-integer (ash integer shift-factor))
569      (setq divisor (ash divisor (- shift-factor)))  ; assume div > num
570      )
571    ;(pprint (list shift-factor scaled-integer divisor))
572    (multiple-value-bind (quotient remainder)(floor scaled-integer divisor)
573      (unless (zerop remainder) ; whats this - tells us there's junk below
574        (setq quotient (logior quotient 1)))
575      ;; why do it return 2 values?
576      (values (float-and-scale-and-round sign quotient (- shift-factor)  short result)))))
577
578
579
580;;; when is (negate-bignum (bignum-ashift-right big)) ; can't negate
581;;; in place cause may get bigger cheaper than (negate-bignum big) - 6
582;;; 0r 8 digits ; 8 longs so win if digits > 7 or negate it on the
583;;; stack
584
585(defun %bignum-dfloat (big &optional result) 
586  (let* ((minusp (bignum-minusp big)))
587    (flet 
588      ((doit (new-big)
589         (let* ((int-len (bignum-integer-length new-big)))
590           (when (>= int-len (- 2047 IEEE-double-float-bias)) ; args?
591             (error (make-condition 'floating-point-overflow 
592                                    :operation 'float :operands (list big))))
593           (if (> int-len 53)
594             (let* ((hi (ldb (byte 25  (- int-len  25)) new-big))
595                    (lo (ldb (byte 28 (- int-len 53)) new-big)))
596               ;(print (list new-big hi lo))
597               (when (and (logbitp (- int-len 54) new-big)  ; round bit
598                          (or (%ilogbitp 0 lo)    ; oddp
599                              ;; or more bits below round
600                              (%i< (one-bignum-factor-of-two new-big) (- int-len 54))))
601                 (if (eq lo #xfffffff)
602                   (setq hi (1+ hi) lo 0)
603                   (setq lo (1+ lo)))
604                 (when (%ilogbitp 25 hi) ; got bigger
605                   (setq int-len (1+ int-len))
606                   (let ((bit (%ilogbitp 0 hi)))
607                     (setq hi (%ilsr 1 hi))
608                     (setq lo (%ilsr 1 lo))
609                     (if bit (setq lo (%ilogior #x8000000 lo))))))
610               (make-float-from-fixnums hi lo (+ IEEE-double-float-bias int-len)(if minusp -1 1) result))
611             (let* ((hi (ldb (byte 25  (- int-len  25)) new-big))
612                    (lobits (min (- int-len 25) 28))
613                    (lo (ldb (byte lobits (- int-len (+ lobits 25))) new-big)))
614               (if (< lobits 28) (setq lo (ash lo (- 28 lobits))))
615               (make-float-from-fixnums hi lo (+ IEEE-double-float-bias int-len) (if minusp -1 1) result))))))
616      (declare (dynamic-extent #'doit))
617      (with-one-negated-bignum-buffer big doit))))
618
619#+32-bit-target
620(defun %bignum-sfloat (big &optional result) 
621  (let* ((minusp (bignum-minusp big)))
622    (flet 
623      ((doit (new-big)
624         (let* ((int-len (bignum-integer-length new-big)))
625           (when (>= int-len (- 255 IEEE-single-float-bias)) ; args?
626             (error (make-condition 'floating-point-overflow 
627                                    :operation 'float :operands (list big 1.0s0))))
628           (if t ;(> int-len IEEE-single-float-digits) ; always true
629             (let* ((lo (ldb (byte IEEE-single-float-digits  (- int-len  IEEE-single-float-digits)) new-big)))
630               (when (and (logbitp (- int-len 25) new-big)  ; round bit
631                          (or (%ilogbitp 0 lo)    ; oddp
632                              ; or more bits below round
633                              (%i< (one-bignum-factor-of-two new-big) (- int-len 25))))
634                 (setq lo (1+ lo))
635                 (when (%ilogbitp 24 lo) ; got bigger
636                   (setq int-len (1+ int-len))
637                   (setq lo (%ilsr 1 lo))))
638               (make-short-float-from-fixnums  lo (+ IEEE-single-float-bias int-len)(if minusp -1 1) result))
639             ))))
640      (declare (dynamic-extent #'doit))
641      (with-one-negated-bignum-buffer big doit))))
642
643
644#+64-bit-target
645(defun %bignum-sfloat (big) 
646  (let* ((minusp (bignum-minusp big)))
647    (flet 
648      ((doit (new-big)
649         (let* ((int-len (bignum-integer-length new-big)))
650           (when (>= int-len (- 255 IEEE-single-float-bias)) ; args?
651             (error (make-condition 'floating-point-overflow 
652                                    :operation 'float :operands (list big 1.0s0))))
653           (if t ;(> int-len IEEE-single-float-digits) ; always true
654             (let* ((lo (ldb (byte IEEE-single-float-digits  (- int-len  IEEE-single-float-digits)) new-big)))
655               (when (and (logbitp (- int-len 25) new-big)  ; round bit
656                          (or (%ilogbitp 0 lo)    ; oddp
657                              ; or more bits below round
658                              (%i< (one-bignum-factor-of-two new-big) (- int-len 25))))
659                 (setq lo (1+ lo))
660                 (when (%ilogbitp 24 lo) ; got bigger
661                   (setq int-len (1+ int-len))
662                   (setq lo (%ilsr 1 lo))))
663               (make-short-float-from-fixnums  lo (+ IEEE-single-float-bias int-len)(if minusp -1 1)))
664             ))))
665      (declare (dynamic-extent #'doit))
666      (with-one-negated-bignum-buffer big doit))))
667
668
669
670
671(defun %fixnum-dfloat (fix &optional result) 
672  (if (eq 0 fix) 
673    (if result (%copy-double-float 0.0d0 result) 0.0d0)
674    (progn
675      (when (not result)(setq result (%make-dfloat)))
676      ; it better return result
677      (%int-to-dfloat fix result))))
678
679
680#+32-bit-target
681(defun %fixnum-sfloat (fix &optional result)
682  (if (eq 0 fix)
683    (if result (%copy-short-float 0.0s0 result) 0.0s0)
684    (%int-to-sfloat! fix (or result (%make-sfloat)))))
685
686#+64-bit-target
687(defun %fixnum-sfloat (fix)
688  (if (eq 0 fix)
689    0.0s0
690    (%int-to-sfloat fix)))
691
692;;; Transcendental functions.
693(defun sin (x)
694  "Return the sine of NUMBER."
695  (if (complexp x)
696    (let* ((r (realpart x))
697           (i (imagpart x)))
698      (complex (* (sin r) (cosh i))
699               (* (cos r) (sinh i))))
700    (if (typep x 'double-float)
701      (%double-float-sin! x (%make-dfloat))
702      #+32-bit-target
703      (target::with-stack-short-floats ((sx x))
704        (%single-float-sin! sx (%make-sfloat)))
705      #+64-bit-target
706      (%single-float-sin (%short-float x)))))
707
708(defun cos (x)
709  "Return the cosine of NUMBER."
710  (if (complexp x)
711    (let* ((r (realpart x))
712           (i (imagpart x)))
713      (complex (* (cos r) (cosh i))
714               (- (* (sin r) (sinh i)))))
715    (if (typep x 'double-float)
716      (%double-float-cos! x (%make-dfloat))
717      #+32-bit-target
718      (target::with-stack-short-floats ((sx x))
719        (%single-float-cos! sx (%make-sfloat)))
720      #+64-bit-target
721      (%single-float-cos (%short-float x)))))
722
723(defun tan (x)
724  "Return the tangent of NUMBER."
725  (if (complexp x)
726    (/ (sin x) (cos x))
727    (if (typep x 'double-float)
728      (%double-float-tan! x (%make-dfloat))
729      #+32-bit-target
730      (target::with-stack-short-floats ((sx x))
731        (%single-float-tan! sx (%make-sfloat)))
732      #+64-bit-target
733      (%single-float-tan (%short-float x))
734      )))
735
736
737;;; Multiply by i (with additional optional scale factor)
738;;; Does the "right thing" with minus zeroes (see CLTL2)
739(defun i* (number &optional (scale 1))
740  (complex (* (- scale) (imagpart number))
741           (* scale (realpart number))))
742
743;;; complex atanh
744(defun %complex-atanh (z)
745  (let* ((rx (realpart z))
746         (ix (imagpart z))
747         (sign (typecase rx
748                 (double-float (%double-float-sign rx))
749                 (short-float (%short-float-sign rx))
750                 (t (minusp rx))))
751         (x rx)
752         (y ix)
753         (y1 (abs y))
754         ra ia)
755    ;; following code requires non-negative x
756    (when sign
757      (setf x (- x))
758      (setf y (- y)))
759    (cond ((> (max x y1) 1.8014399e+16)
760           ;; large value escape
761           (setq ra (if (> x y1)
762                      (let ((r (/ y x)))
763                        (/ (/ x) (1+ (* r r))))
764                      (let ((r (/ x y)))
765                        (/ (/ r y) (1+ (* r r))))))
766           (setq ia (typecase y
767                      (double-float (float-sign y double-float-half-pi))
768                      (single-float (float-sign y single-float-half-pi))
769                      (t (if (minusp y) #.(- single-float-half-pi) single-float-half-pi)))))
770          ((= x 1)
771           (setq ra (if (< y1 1e-9)
772                      (/ (log-e (/ 2 y1)) 2)
773                      (/ (log1+ (/ 4 (* y y))) 4)))
774           (setq ia (/ (atan (/ 2 y) -1) 2)))
775          (t
776           (let ((r2 (+ (* x x) (* y y))))
777             (setq ra (/ (log1+ (/ (* -4 x) (1+ (+ (* 2 x) r2)))) -4))
778             (setq ia (/ (atan (* 2 y) (- 1 r2)) 2)))))
779    ;; fixup signs, with special case for real arguments
780    (cond (sign
781           (setq ra (- ra))
782           (when (typep z 'complex)
783             (setq ia (- ia))))
784          (t
785           (unless (typep z 'complex)
786             (setq ia (- ia)))))
787    (complex ra ia)))
788
789(defun atan (y &optional (x nil x-p))
790  "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
791  (cond (x-p
792         (cond ((or (typep x 'double-float)
793                    (typep y 'double-float))
794                (with-stack-double-floats ((dy y)
795                                           (dx x))
796                  (%df-atan2 dy dx)))
797               (t
798                (when (and (rationalp x) (rationalp y))
799                  ;; rescale arguments so that the maximum absolute value is 1
800                  (let ((x1 (abs x)) (y1 (abs y)))
801                    (cond ((> y1 x1)
802                           (setf x (/ x y1))
803                           (setf y (signum y)))
804                          ((not (zerop x))
805                           (setf y (/ y x1))
806                           (setf x (signum x))))))
807                #+32-bit-target
808                (target::with-stack-short-floats ((sy y)
809                                                  (sx x))
810                  (%sf-atan2! sy sx))
811                #+64-bit-target
812                (%sf-atan2 (%short-float y) (%short-float x)))))
813        ((typep y 'double-float)
814         (%double-float-atan! y (%make-dfloat)))
815        ((typep y 'single-float)
816         #+32-bit-target
817         (%single-float-atan! y (%make-sfloat))
818         #+64-bit-target
819         (%single-float-atan y))
820        ((typep y 'rational)
821         (cond ((<= (abs y) most-positive-short-float)
822                #+32-bit-target
823                (target::with-stack-short-floats ((sy y))
824                  (%single-float-atan! sy (%make-sfloat)))
825                #+64-bit-target
826                (%single-float-atan (%short-float y)))
827               ((minusp y)
828                #.(- single-float-half-pi))
829               (t
830                single-float-half-pi)))
831        (t
832         (let ((r (realpart y))
833               (i (imagpart y)))
834           (if (zerop i)
835             (complex (atan r) i)
836             (i* (%complex-atanh (complex (- i) r)) -1))))))
837
838
839
840(defun log (x &optional (b nil b-p))
841  "Return the logarithm of NUMBER in the base BASE, which defaults to e."
842  (if b-p
843    (if (zerop b)
844      (if (zerop x)
845        (report-bad-arg x '(not (satisfies zerop) ))
846        ;; ** CORRECT THE CONTAGION for complex args **
847        (+ 0 (* x b)))
848      ;; do the float/rational contagion before the division
849      ;; but take care with negative zeroes
850      (let ((x1 (realpart x))
851            (b1 (realpart b)))
852        (if (and (typep x1 'float)
853                 (typep b1 'float))
854          (/ (log-e (* x (float 1.0 b1)))
855             (log-e (* b (float 1.0 x1))))
856          (let ((r (/ (cond ((typep x 'rational)
857                             (%rational-log x 1.0d0))
858                            ((typep x1 'rational)
859                             (%rational-complex-log x 1.0d0))
860                            (t
861                             (log-e (* x 1.0d0))))
862                      (cond ((typep b 'rational)
863                             (%rational-log b 1.0d0))
864                            ((typep b1 'rational)
865                             (%rational-complex-log b 1.0d0))
866                            (t
867                             (log-e (* b 1.0d0)))))))
868            (cond ((or (typep x1 'double-float)
869                       (typep b1 'double-float))
870                   r)
871                  ((complexp r)
872                   (complex (%short-float (realpart r))
873                            (%short-float (imagpart r))))
874                  (t
875                   (%short-float r)))))))
876    (log-e x)))
877
878
879
880(defun log-e (x)
881   (cond
882     ((typep x 'double-float)
883      (if (%double-float-sign x)
884        (with-stack-double-floats ((dx x))
885          (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi))
886        (%double-float-log! x (%make-dfloat))))
887    ((typep x 'short-float)
888     #+32-bit-target
889     (if (%short-float-sign x)
890       (target::with-stack-short-floats ((sx x))
891         (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat))
892                  single-float-pi))
893       (%single-float-log! x (%make-sfloat)))
894     #+64-bit-target
895     (if (%short-float-sign x)
896       (complex (%single-float-log (%short-float-abs (%short-float x)))
897                single-float-pi)
898       (%single-float-log (%short-float x))))
899    ((typep x 'complex)
900     (if (typep (realpart x) 'rational)
901       (%rational-complex-log x 1.0s0)
902       ;; take care that intermediate results do not overflow/underflow:
903       ;; pre-scale argument by an appropriate power of two;
904       ;; we only need to scale for very large and very small values -
905       ;;  hence the various 'magic' numbers (values may not be too
906       ;;  critical but do depend on the sqrt update to fix abs's operation)
907       (let ((m (max (abs (realpart x))
908                     (abs (imagpart x))))
909             (log-s 0)
910             (s 1))
911         (if (typep m 'short-float)
912           (let ((expon (- (%short-float-exp m) IEEE-single-float-bias)))
913             (cond ((> expon 126)
914                    (setq log-s double-float-log2^23)
915                    (setq s #.(ash 1 23)))
916                   ((< expon -124)
917                    (setq log-s #.(- double-float-log2^23))
918                    (setq s #.(/ 1.0s0 (ash 1 23))))))
919           (let ((expon (- (%double-float-exp m) IEEE-double-float-bias)))
920             (cond ((> expon 1022)
921                    (setq log-s double-float-log2^23)
922                    (setq s #.(ash 1 23)))
923                   ((< expon -1020)
924                    (setq log-s #.(- double-float-log2^23))
925                    (setq s #.(/ 1.0d0 (ash 1 23)))))))
926         (if (eql s 1)
927           (complex (log-abs x) (phase x))
928           (let ((temp (log-abs (/ x s))))
929             (complex (float (+ log-s temp) temp) (phase x)))))))
930    (t
931     (%rational-log x 1.0s0))))
932
933;;; helpers for rational log
934(defun %rational-log (x prototype)
935  (cond ((minusp x)
936         (complex (%rational-log (- x) prototype)
937                  (if (typep prototype 'short-float)
938                    single-float-pi
939                    pi)))
940        ((bignump x)
941         ;(let* ((base1 3)
942         ;       (guess (floor (1- (integer-length x))
943         ;                     (log base1 2)))
944         ;       (guess1 (* guess (log-e base1))))
945         ;  (+ guess1 (log-e (/ x (expt base1 guess)))))
946         ; Using base1 = 2 is *much* faster. Was there a reason for 3?
947         (let* ((guess (1- (integer-length x)))
948                (guess1 (* guess double-float-log2)))
949           (float (+ guess1 (log-e (float (/ x (ash 1 guess)) 1.0d0))) prototype)))
950        ((and (ratiop x)
951              ;; Rational arguments near +1 can be specified with great precision: don't lose it
952              (cond ((< 0.5 x 3)
953                     (log1+ (float (- x 1) prototype)))
954                    (
955                     ;; Escape out small values as well as large
956                     (or (> x most-positive-short-float)
957                         (< x least-positive-normalized-short-float))
958                     ;; avoid loss of precision due to subtracting logs of numerator and denominator
959                     (let* ((n (%numerator x))
960                            (d (%denominator x))
961                            (sn (1- (integer-length n)))
962                            (sd (1- (integer-length d))))
963                       (float (+ (* (- sn sd) double-float-log2)
964                                 (- (log1+ (float (1- (/ n (ash 1 sn))) 1.0d0))
965                                    (log1+ (float (1- (/ d (ash 1 sd))) 1.0d0))))
966                              prototype))))))
967        ((typep prototype 'short-float)
968         #+32-bit-target
969         (target::with-stack-short-floats ((sx x))
970           (%single-float-log! sx (%make-sfloat)))
971         #+64-bit-target
972         (%single-float-log (%short-float x)))
973        (t
974         (with-stack-double-floats ((dx x))
975           (%double-float-log! dx (%make-dfloat))))))
976
977;;; (log1+ x) = (log (1+ x))
978;;; but is much more precise near x = 0
979(defun log1+ (x)
980  ;;(cond ((typep x 'complex)
981  ;;      (let ((r (realpart x))
982  ;;            (i (imagpart x)))
983  ;;        (if (and (< (abs r) 0.5)
984  ;;                 (< (abs i) 3))
985  ;;          (let* ((n (+ (* r (+ 2 r)) (* i i)))
986  ;;                 (d (1+ (sqrt (1+ n)))))
987  ;;            (complex (log1+ (/ n d)) (atan i (1+ r))))
988  ;;         (log (1+ x)))))
989  ;;     (t
990  (if (and (typep x 'ratio)
991           (< -0.5 x 2))
992    (setq x (%short-float x)))
993  (let ((y (1+ x)))
994    (if (eql y x)
995      (log-e y)
996      (let ((e (1- y)))
997        (if (zerop e)
998          (* x 1.0)
999          (- (log-e y) (/ (- e x) y)))))))
1000
1001;;; helper for complex log
1002;;; uses more careful approach when (abs x) is near 1
1003(defun log-abs (x)
1004  (let ((a (abs x)))
1005    (if (< 0.5 a 3)
1006      (let* ((r (realpart x))
1007             (i (imagpart x))
1008             (n (if (> (abs r) (abs i))
1009                  (+ (* (1+ r) (1- r)) (* i i))
1010                  (+ (* r r) (* (1+ i) (1- i))))))
1011        (log1+ (/ n (1+ a))))
1012      (log-e a))))
1013
1014(defun %rational-complex-log (x prototype &aux ra ia)
1015  (let* ((rx (realpart x))
1016         (ix (imagpart x))
1017         (x (abs rx))
1018         (y (abs ix)))
1019    (if (> y x)
1020      (let ((r (float (/ rx y) 1.0d0)))
1021        (setq ra (+ (%rational-log y 1.0d0)
1022                    (/ (log1+ (* r r)) 2)))
1023        (setq ia (atan (if (minusp ix) -1.0d0 1.0d0) r)))
1024      (let ((r (float (/ ix x) 1.0d0)))
1025        (setq ra (+ (%rational-log x 1.0d0)
1026                    (/ (log1+ (* r r)) 2)))
1027        (setq ia (atan r (if (minusp rx) -1.0d0 1.0d0)))))
1028    (if (typep prototype 'short-float)
1029      (complex (%short-float ra) (%short-float ia))
1030      (complex ra ia))))
1031
1032(defun exp (x)
1033  "Return e raised to the power NUMBER."
1034  (typecase x
1035    (complex (* (exp (realpart x)) (cis (imagpart x))))
1036    (double-float (%double-float-exp! x (%make-dfloat)))
1037    (t
1038     #+32-bit-target
1039     (target::with-stack-short-floats ((sx x))
1040       (%single-float-exp! sx (%make-sfloat)))
1041     #+64-bit-target
1042     (%single-float-exp (%short-float x)))))
1043
1044
1045(defun positive-realpart-p (n)
1046  (> (realpart n) 0))
1047
1048(defun expt (b e)
1049  "Return BASE raised to the POWER."
1050  (cond ((zerop e) (1+ (* b e)))
1051        ((integerp e)
1052         (if (minusp e) (/ 1 (%integer-power b (- e))) (%integer-power b e)))
1053        ((zerop b)
1054         (if (plusp (realpart e)) b (report-bad-arg e '(satisfies positive-realpart-p))))
1055        ((and (realp b) (plusp b) (realp e))
1056         (if (or (typep b 'double-float)
1057                 (typep e 'double-float))
1058           (with-stack-double-floats ((b1 b)
1059                                      (e1 e))
1060             (%double-float-expt! b1 e1 (%make-dfloat)))
1061           #+32-bit-target
1062           (target::with-stack-short-floats ((b1 b)
1063                                     (e1 e))
1064             (%single-float-expt! b1 e1 (%make-sfloat)))
1065           #+64-bit-target
1066           (%single-float-expt (%short-float b) (%short-float e))
1067           ))
1068        ((typep (realpart e) 'double-float)
1069         ;; Avoid intermediate single-float result from LOG
1070         (let ((promoted-base (* 1d0 b)))
1071           (exp (* e (log promoted-base)))))
1072        (t (exp (* e (log b))))))
1073
1074
1075       
1076(defun sqrt (x &aux a b)
1077  "Return the square root of NUMBER."
1078  (cond ((zerop x) x)
1079        ((complexp x)
1080         (let ((rx (realpart x))
1081               (ix (imagpart x)))
1082           (cond ((rationalp rx)
1083                  (if (zerop rx)
1084                    (let ((s (sqrt (/ (abs ix) 2))))
1085                      (complex s (if (minusp ix) (- s) s)))
1086                    (let* ((s (+ (* rx rx) (* ix ix)))
1087                           (d (if (ratiop s)
1088                                (/ (isqrt (%numerator s))
1089                                   (isqrt (%denominator s)))
1090                                (isqrt s))))
1091                      (unless (eql s (* d d))
1092                        (setf d (%double-float-hypot (%double-float rx)
1093                                                     (%double-float ix))))
1094                      (cond ((minusp rx)
1095                             (setq b (sqrt (/ (- d rx) 2)))
1096                             (when (minusp ix)
1097                               (setq b (- b)))
1098                             (setq a (/ ix (+ b b))))
1099                            (t
1100                             (setq a (sqrt (/ (+ rx d) 2)))
1101                             (setq b (/ ix (+ a a)))))
1102                      (if (rationalp a)
1103                        (complex a b)
1104                        (complex (%short-float a) (%short-float b))))))
1105                 ((minusp rx)
1106                  (if (zerop ix)
1107                    (complex 0 (float-sign ix (sqrt (- rx))))
1108                    (let ((shift (cond ((< rx -1) -3)
1109                                       ((and (> rx -5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25)
1110                                       (t -1))))
1111                      (setq rx (scale-float rx shift))
1112                      (let ((s (fsqrt (- (abs (complex rx (scale-float ix shift))) rx))))
1113                        (setq b (scale-float s (ash (- -1 shift) -1)))
1114                        (when (minusp ix)
1115                          (setq b (- b)))
1116                        (setq a (/ ix (scale-float b 1)))
1117                        (complex a b)))))
1118                 (t
1119                  (if (zerop ix)
1120                    (complex (sqrt rx) ix)
1121                    (let ((shift (cond ((> rx 1) -3)
1122                                       ((and (< rx 5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25)
1123                                       (t -1))))
1124                      (setq rx (scale-float rx shift))
1125                      (let ((s (fsqrt (+ rx (abs (complex rx (scale-float ix shift)))))))
1126                        (setq a (scale-float s (ash (- -1 shift) -1)))
1127                        (setq b (/ ix (scale-float a 1)))
1128                        (complex a b))))))))
1129        ((minusp x) (complex 0 (sqrt (- x))))
1130        ((floatp x)
1131         (fsqrt x))
1132        ((and (integerp x) (eql x (* (setq a (isqrt x)) a))) a)
1133        ((and (ratiop x)
1134              (let ((n (numerator x))
1135                    d)
1136                (and (eql n (* (setq a (isqrt n)) a))
1137                     (eql (setq d (denominator x))
1138                          (* (setq b (isqrt d)) b)))))
1139         (/ a b))         
1140        (t
1141         (float (fsqrt (float x 0.0d0)) 1.0s0))))
1142
1143
1144
1145(defun asin (x)
1146  "Return the arc sine of NUMBER."
1147  (number-case x
1148    (complex
1149      (let ((sqrt-1-x (sqrt (- 1 x)))
1150            (sqrt-1+x (sqrt (+ 1 x))))
1151        (complex (atan (realpart x)
1152                       (realpart (* sqrt-1-x sqrt-1+x)))
1153                 (asinh (imagpart (* (conjugate sqrt-1-x)
1154                                     sqrt-1+x))))))
1155    (double-float
1156     (locally (declare (type double-float x))
1157       (if (and (<= -1.0d0 x)
1158                (<= x 1.0d0))
1159         (%double-float-asin! x (%make-dfloat))
1160         (let* ((temp (+ (complex -0.0d0 x)
1161                         (sqrt (- 1.0d0 (the double-float (* x x)))))))
1162           (complex (phase temp) (- (log (abs temp))))))))
1163    ((short-float rational)
1164     #+32-bit-target
1165     (let* ((x1 (%make-sfloat)))
1166       (declare (dynamic-extent x1))
1167       (if (and (realp x) 
1168                (<= -1.0s0 (setq x (%short-float x x1)))
1169                (<= x 1.0s0))
1170         (%single-float-asin! x1 (%make-sfloat))
1171         (progn
1172           (setq x (+ (complex (- (imagpart x)) (realpart x))
1173                      (sqrt (- 1 (* x x)))))
1174           (complex (phase x) (- (log (abs x)))))))
1175     #+64-bit-target
1176     (if (and (realp x) 
1177              (<= -1.0s0 (setq x (%short-float x)))
1178              (<= x 1.0s0))
1179         (%single-float-asin x)
1180         (progn
1181           (setq x (+ (complex (- (imagpart x)) (realpart x))
1182                      (sqrt (- 1 (* x x)))))
1183           (complex (phase x) (- (log (abs x))))))
1184     )))
1185
1186
1187(defun acos (x)
1188  "Return the arc cosine of NUMBER."
1189  (number-case x
1190    (complex
1191     (if (typep (realpart x) 'double-float)
1192       (- double-float-half-pi (asin x))
1193       (- single-float-half-pi (asin x))))
1194    (double-float
1195     (locally (declare (type double-float x))
1196       (if (and (<= -1.0d0 x)
1197                (<= x 1.0d0))
1198         (%double-float-acos! x (%make-dfloat))
1199         (- double-float-half-pi (asin x)))))
1200    ((short-float rational)
1201     #+32-bit-target
1202     (target::with-stack-short-floats ((sx x))
1203        (locally
1204            (declare (type short-float sx))
1205          (if (and (<= -1.0s0 sx)
1206                   (<= sx 1.0s0))
1207            (%single-float-acos! sx (%make-sfloat))
1208            (- single-float-half-pi (asin sx)))))
1209     #+64-bit-target
1210     (let* ((sx (%short-float x)))
1211       (declare (type short-float sx))
1212       (if (and (<= -1.0s0 sx)
1213                (<= sx 1.0s0))
1214         (%single-float-acos sx)
1215         (- single-float-half-pi (asin sx))))
1216     )))
1217
1218
1219(defun fsqrt (x)
1220  (etypecase x
1221    (double-float (%double-float-sqrt! x (%make-dfloat)))
1222    (single-float
1223     #+32-bit-target
1224     (%single-float-sqrt! x (%make-sfloat))
1225     #+64-bit-target
1226     (%single-float-sqrt x))))
1227
1228
1229
1230(defun %df-atan2 (y x)
1231  (if (zerop x)
1232    (if (zerop y)
1233      (if (plusp (float-sign x))
1234        (if (eql y -0.0d0)
1235          -0.0d0
1236          0.0d0)
1237        (float-sign y pi))
1238      (float-sign y double-float-half-pi))
1239    (%double-float-atan2! y x (%make-dfloat))))
1240
1241#+32-bit-target
1242(defun %sf-atan2! (y x)
1243  (if (zerop x)
1244    (if (zerop y)
1245      (if (plusp (float-sign x))
1246        ;; Don't return Y (which may be stack-consed) here.
1247        ;; We know that (ZEROP Y) is true, so:
1248        (if (eql y -0.0s0)
1249          -0.0s0
1250          0.0s0)
1251        (float-sign y single-float-pi))
1252      (float-sign y single-float-half-pi))
1253    (%single-float-atan2! y x (%make-sfloat))))
1254
1255#+64-bit-target
1256(defun %sf-atan2 (y x)
1257  (if (zerop x)
1258    (if (zerop y)
1259      (if (plusp (float-sign x))
1260        y
1261        (float-sign y single-float-pi))
1262      (float-sign y single-float-half-pi))
1263    (%single-float-atan2 y x)))
1264
1265#+64-bit-target
1266(defun %short-float-exp (n)
1267  (let* ((bits (single-float-bits n)))
1268    (declare (type (unsigned-byte 32) bits))
1269    (ldb (byte IEEE-single-float-exponent-width IEEE-single-float-exponent-offset) bits)))
1270
1271
1272#+64-bit-target
1273(defun set-%short-float-exp (float exp)
1274  (host-single-float-from-unsigned-byte-32
1275   (dpb exp
1276        (byte IEEE-single-float-exponent-width
1277              IEEE-single-float-exponent-offset)
1278        (the (unsigned-byte 32) (single-float-bits float)))))
1279
1280#+64-bit-target
1281(defun %%scale-sfloat (float int)
1282  (* (the single-float float)
1283     (the single-float (host-single-float-from-unsigned-byte-32
1284                        (dpb int
1285                             (byte IEEE-single-float-exponent-width
1286                                   IEEE-single-float-exponent-offset)
1287                             0)))))
1288
1289#+64-bit-target
1290(defun %double-float-exp (n)
1291  (let* ((highword (double-float-bits n)))
1292    (declare (fixnum highword))
1293    (logand (1- (ash 1 IEEE-double-float-exponent-width))
1294            (ash highword (- (- IEEE-double-float-exponent-offset 32))))))
1295
1296#+64-bit-target
1297(defun set-%double-float-exp (float exp)
1298  (let* ((highword (double-float-bits float)))
1299    (declare (fixnum highword))
1300    (setf (uvref float target::double-float.val-high-cell)
1301          (dpb exp
1302               (byte IEEE-double-float-exponent-width
1303                     (- IEEE-double-float-exponent-offset 32))
1304               highword))
1305    exp))
1306
1307#+64-bit-target
1308(defun %integer-decode-double-float (f)
1309  (multiple-value-bind (hiword loword) (double-float-bits f)
1310    (declare (type (unsigned-byte 32) hiword loword))
1311    (let* ((exp (ldb (byte IEEE-double-float-exponent-width
1312                           (- IEEE-double-float-exponent-offset 32))
1313                     hiword))
1314           (mantissa (logior
1315                      (the fixnum
1316                        (dpb (ldb (byte (- IEEE-double-float-mantissa-width 32)
1317                                        IEEE-double-float-mantissa-offset)
1318                                  hiword)
1319                             (byte (- IEEE-double-float-mantissa-width 32)
1320                                   32)
1321                             loword))
1322                      (if (zerop exp)
1323                        0
1324                        (ash 1 IEEE-double-float-hidden-bit))))
1325           (sign (if (logbitp 31 hiword) -1 1)))
1326      (declare (fixnum exp mantissa sign))
1327      (values (ldb (byte 25 28) mantissa)
1328              (ldb (byte 28 0) mantissa)
1329              exp
1330              sign))))
1331
1332;;; end of l0-float.lisp
Note: See TracBrowser for help on using the repository browser.