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

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

Propagate recent trunk changes.

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