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

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

Implement the rest of Dave Findlay's "reference" math library.
Fixes ticket:869 in the trunk.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 63.9 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                    (n (if (> (abs i) 20)
748                         (* 4 (exp (* -2 (abs i))))
749                         (let ((c (cosh i)))
750                           (/ (* c c))))))
751               (complex (/ (* n tx) d)
752                        (/ (* ty (1+ tx2)) d))))))
753        ((or (typep x 'ratio)
754             (> (abs x) two^23))
755         (let ((c (%extended-cis x)))
756           (if (typep x 'double-float)
757             (/ (imagpart c) (realpart c))
758             (%short-float (/ (imagpart c) (realpart c))))))
759        ((typep x 'double-float)
760         (%double-float-tan! x (%make-dfloat)))
761        (t
762         #+32-bit-target
763         (target::with-stack-short-floats ((sx x))
764           (%single-float-tan! sx (%make-sfloat)))
765         #+64-bit-target
766         (%single-float-tan (%short-float x))
767         )))
768
769
770;;; Helper function for sin/cos/tan for ratio or large arguments
771;;; (Will become inaccurate for ridiculously large arguments though)
772;;; Does not assume that float library returns accurate values for large arguments
773;;; (many don't)
774
775;;; hexdecimal representations of pi at various precisions
776(defconstant pi-vector
777  #(#x3243F6A8885A308D313198A2E0
778    #x3243F6A8885A308D313198A2E03707344A4093822299F31D008
779    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D
780    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC
781    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B5470
782    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310B
783    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB
784    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045
785    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70
786    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871
787    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D9
788    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD588215
789    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D
790    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085
791    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB
792    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3E
793    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3ED71577C1BD314B2778AF2FDA5
794    ))
795
796(defun %extended-cis (x)
797  (let (size pi-size)
798    (typecase x
799      (integer (setq size (1- (integer-length (abs x)))))
800      (ratio (setq size (- (integer-length (abs (numerator x)))
801                           (integer-length (denominator x)))))
802      (float (multiple-value-bind (mantissa exponent sign)
803                                  (integer-decode-float x)
804               (setq size (+ (1- (integer-length mantissa)) exponent))
805               (setq x (* sign mantissa (expt 2 exponent))))))
806    (setq pi-size (ceiling (+ size 64) 100))
807    (cond ((< pi-size 1) (setq pi-size 1))
808          ((> pi-size 17) (setq pi-size 17)))
809    (let* ((2pi-approx (/ (aref pi-vector (1- pi-size))
810                          (ash 1 (1- (* 100 pi-size)))))
811           (reduced-x (rem x 2pi-approx))
812           (x0 (float reduced-x 1.0d0))
813           (x1 (- reduced-x (rational x0))))
814      (* (cis x0) (cis (float x1 1.0d0))))))
815
816
817;;; Multiply by i (with additional optional scale factor)
818;;; Does the "right thing" with minus zeroes (see CLTL2)
819(defun i* (number &optional (scale 1))
820  (complex (* (- scale) (imagpart number))
821           (* scale (realpart number))))
822
823;;; complex atanh
824(defun %complex-atanh (z)
825  (let* ((rx (realpart z))
826         (ix (imagpart z))
827         (sign (typecase rx
828                 (double-float (%double-float-sign rx))
829                 (short-float (%short-float-sign rx))
830                 (t (minusp rx))))
831         (x rx)
832         (y ix)
833         (y1 (abs y))
834         ra ia)
835    ;; following code requires non-negative x
836    (when sign
837      (setf x (- x))
838      (setf y (- y)))
839    (cond ((> (max x y1) 1.8014399e+16)
840           ;; large value escape
841           (setq ra (if (> x y1)
842                      (let ((r (/ y x)))
843                        (/ (/ x) (1+ (* r r))))
844                      (let ((r (/ x y)))
845                        (/ (/ r y) (1+ (* r r))))))
846           (setq ia (typecase y
847                      (double-float (float-sign y double-float-half-pi))
848                      (single-float (float-sign y single-float-half-pi))
849                      (t (if (minusp y) #.(- single-float-half-pi) single-float-half-pi)))))
850          ((= x 1)
851           (setq ra (if (< y1 1e-9)
852                      (/ (log-e (/ 2 y1)) 2)
853                      (/ (log1+ (/ 4 (* y y))) 4)))
854           (setq ia (/ (atan (/ 2 y) -1) 2)))
855          (t
856           (let ((r2 (+ (* x x) (* y y))))
857             (setq ra (/ (log1+ (/ (* -4 x) (1+ (+ (* 2 x) r2)))) -4))
858             (setq ia (/ (atan (* 2 y) (- 1 r2)) 2)))))
859    ;; fixup signs, with special case for real arguments
860    (cond (sign
861           (setq ra (- ra))
862           (when (typep z 'complex)
863             (setq ia (- ia))))
864          (t
865           (unless (typep z 'complex)
866             (setq ia (- ia)))))
867    (complex ra ia)))
868
869(defun atan (y &optional (x nil x-p))
870  "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
871  (cond (x-p
872         (cond ((or (typep x 'double-float)
873                    (typep y 'double-float))
874                (with-stack-double-floats ((dy y)
875                                           (dx x))
876                  (%df-atan2 dy dx)))
877               (t
878                (when (and (rationalp x) (rationalp y))
879                  ;; rescale arguments so that the maximum absolute value is 1
880                  (let ((x1 (abs x)) (y1 (abs y)))
881                    (cond ((> y1 x1)
882                           (setf x (/ x y1))
883                           (setf y (signum y)))
884                          ((not (zerop x))
885                           (setf y (/ y x1))
886                           (setf x (signum x))))))
887                #+32-bit-target
888                (target::with-stack-short-floats ((sy y)
889                                                  (sx x))
890                  (%sf-atan2! sy sx))
891                #+64-bit-target
892                (%sf-atan2 (%short-float y) (%short-float x)))))
893        ((typep y 'double-float)
894         (%double-float-atan! y (%make-dfloat)))
895        ((typep y 'single-float)
896         #+32-bit-target
897         (%single-float-atan! y (%make-sfloat))
898         #+64-bit-target
899         (%single-float-atan y))
900        ((typep y 'rational)
901         (cond ((<= (abs y) most-positive-short-float)
902                #+32-bit-target
903                (target::with-stack-short-floats ((sy y))
904                  (%single-float-atan! sy (%make-sfloat)))
905                #+64-bit-target
906                (%single-float-atan (%short-float y)))
907               ((minusp y)
908                #.(- single-float-half-pi))
909               (t
910                single-float-half-pi)))
911        (t
912         (let ((r (realpart y))
913               (i (imagpart y)))
914           (if (zerop i)
915             (complex (atan r) i)
916             (i* (%complex-atanh (complex (- i) r)) -1))))))
917
918
919
920(defun log (x &optional (b nil b-p))
921  "Return the logarithm of NUMBER in the base BASE, which defaults to e."
922  (if b-p
923    (if (zerop b)
924      (if (zerop x)
925        (report-bad-arg x '(not (satisfies zerop) ))
926        ;; ** CORRECT THE CONTAGION for complex args **
927        (+ 0 (* x b)))
928      ;; do the float/rational contagion before the division
929      ;; but take care with negative zeroes
930      (let ((x1 (realpart x))
931            (b1 (realpart b)))
932        (if (and (typep x1 'float)
933                 (typep b1 'float))
934          (/ (log-e (* x (float 1.0 b1)))
935             (log-e (* b (float 1.0 x1))))
936          (let ((r (/ (cond ((typep x 'rational)
937                             (%rational-log x 1.0d0))
938                            ((typep x1 'rational)
939                             (%rational-complex-log x 1.0d0))
940                            (t
941                             (log-e (* x 1.0d0))))
942                      (cond ((typep b 'rational)
943                             (%rational-log b 1.0d0))
944                            ((typep b1 'rational)
945                             (%rational-complex-log b 1.0d0))
946                            (t
947                             (log-e (* b 1.0d0)))))))
948            (cond ((or (typep x1 'double-float)
949                       (typep b1 'double-float))
950                   r)
951                  ((complexp r)
952                   (complex (%short-float (realpart r))
953                            (%short-float (imagpart r))))
954                  (t
955                   (%short-float r)))))))
956    (log-e x)))
957
958
959
960(defun log-e (x)
961   (cond
962     ((typep x 'double-float)
963      (if (%double-float-sign x)
964        (with-stack-double-floats ((dx x))
965          (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi))
966        (%double-float-log! x (%make-dfloat))))
967    ((typep x 'short-float)
968     #+32-bit-target
969     (if (%short-float-sign x)
970       (target::with-stack-short-floats ((sx x))
971         (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat))
972                  single-float-pi))
973       (%single-float-log! x (%make-sfloat)))
974     #+64-bit-target
975     (if (%short-float-sign x)
976       (complex (%single-float-log (%short-float-abs (%short-float x)))
977                single-float-pi)
978       (%single-float-log (%short-float x))))
979    ((typep x 'complex)
980     (if (typep (realpart x) 'rational)
981       (%rational-complex-log x 1.0s0)
982       ;; take care that intermediate results do not overflow/underflow:
983       ;; pre-scale argument by an appropriate power of two;
984       ;; we only need to scale for very large and very small values -
985       ;;  hence the various 'magic' numbers (values may not be too
986       ;;  critical but do depend on the sqrt update to fix abs's operation)
987       (let ((m (max (abs (realpart x))
988                     (abs (imagpart x))))
989             (log-s 0)
990             (s 1))
991         (if (typep m 'short-float)
992           (let ((expon (- (%short-float-exp m) IEEE-single-float-bias)))
993             (cond ((> expon 126)
994                    (setq log-s double-float-log2^23)
995                    (setq s #.(ash 1 23)))
996                   ((< expon -124)
997                    (setq log-s #.(- double-float-log2^23))
998                    (setq s #.(/ 1.0s0 (ash 1 23))))))
999           (let ((expon (- (%double-float-exp m) IEEE-double-float-bias)))
1000             (cond ((> expon 1022)
1001                    (setq log-s double-float-log2^23)
1002                    (setq s #.(ash 1 23)))
1003                   ((< expon -1020)
1004                    (setq log-s #.(- double-float-log2^23))
1005                    (setq s #.(/ 1.0d0 (ash 1 23)))))))
1006         (if (eql s 1)
1007           (complex (log-abs x) (phase x))
1008           (let ((temp (log-abs (/ x s))))
1009             (complex (float (+ log-s temp) temp) (phase x)))))))
1010    (t
1011     (%rational-log x 1.0s0))))
1012
1013;;; helpers for rational log
1014(defun %rational-log (x prototype)
1015  (cond ((minusp x)
1016         (complex (%rational-log (- x) prototype)
1017                  (if (typep prototype 'short-float)
1018                    single-float-pi
1019                    pi)))
1020        ((bignump x)
1021         ;(let* ((base1 3)
1022         ;       (guess (floor (1- (integer-length x))
1023         ;                     (log base1 2)))
1024         ;       (guess1 (* guess (log-e base1))))
1025         ;  (+ guess1 (log-e (/ x (expt base1 guess)))))
1026         ; Using base1 = 2 is *much* faster. Was there a reason for 3?
1027         (let* ((guess (1- (integer-length x)))
1028                (guess1 (* guess double-float-log2)))
1029           (float (+ guess1 (log-e (float (/ x (ash 1 guess)) 1.0d0))) prototype)))
1030        ((and (ratiop x)
1031              ;; Rational arguments near +1 can be specified with great precision: don't lose it
1032              (cond ((< 0.5 x 3)
1033                     (log1+ (float (- x 1) prototype)))
1034                    (
1035                     ;; Escape out small values as well as large
1036                     (or (> x most-positive-short-float)
1037                         (< x least-positive-normalized-short-float))
1038                     ;; avoid loss of precision due to subtracting logs of numerator and denominator
1039                     (let* ((n (%numerator x))
1040                            (d (%denominator x))
1041                            (sn (1- (integer-length n)))
1042                            (sd (1- (integer-length d))))
1043                       (float (+ (* (- sn sd) double-float-log2)
1044                                 (- (log1+ (float (1- (/ n (ash 1 sn))) 1.0d0))
1045                                    (log1+ (float (1- (/ d (ash 1 sd))) 1.0d0))))
1046                              prototype))))))
1047        ((typep prototype 'short-float)
1048         #+32-bit-target
1049         (target::with-stack-short-floats ((sx x))
1050           (%single-float-log! sx (%make-sfloat)))
1051         #+64-bit-target
1052         (%single-float-log (%short-float x)))
1053        (t
1054         (with-stack-double-floats ((dx x))
1055           (%double-float-log! dx (%make-dfloat))))))
1056
1057;;; (log1+ x) = (log (1+ x))
1058;;; but is much more precise near x = 0
1059(defun log1+ (x)
1060  ;;(cond ((typep x 'complex)
1061  ;;      (let ((r (realpart x))
1062  ;;            (i (imagpart x)))
1063  ;;        (if (and (< (abs r) 0.5)
1064  ;;                 (< (abs i) 3))
1065  ;;          (let* ((n (+ (* r (+ 2 r)) (* i i)))
1066  ;;                 (d (1+ (sqrt (1+ n)))))
1067  ;;            (complex (log1+ (/ n d)) (atan i (1+ r))))
1068  ;;         (log (1+ x)))))
1069  ;;     (t
1070  (if (and (typep x 'ratio)
1071           (< -0.5 x 2))
1072    (setq x (%short-float x)))
1073  (let ((y (1+ x)))
1074    (if (eql y x)
1075      (log-e y)
1076      (let ((e (1- y)))
1077        (if (zerop e)
1078          (* x 1.0)
1079          (- (log-e y) (/ (- e x) y)))))))
1080
1081;;; helper for complex log
1082;;; uses more careful approach when (abs x) is near 1
1083(defun log-abs (x)
1084  (let ((a (abs x)))
1085    (if (< 0.5 a 3)
1086      (let* ((r (realpart x))
1087             (i (imagpart x))
1088             (n (if (> (abs r) (abs i))
1089                  (+ (* (1+ r) (1- r)) (* i i))
1090                  (+ (* r r) (* (1+ i) (1- i))))))
1091        (log1+ (/ n (1+ a))))
1092      (log-e a))))
1093
1094(defun %rational-complex-log (x prototype &aux ra ia)
1095  (let* ((rx (realpart x))
1096         (ix (imagpart x))
1097         (x (abs rx))
1098         (y (abs ix)))
1099    (if (> y x)
1100      (let ((r (float (/ rx y) 1.0d0)))
1101        (setq ra (+ (%rational-log y 1.0d0)
1102                    (/ (log1+ (* r r)) 2)))
1103        (setq ia (atan (if (minusp ix) -1.0d0 1.0d0) r)))
1104      (let ((r (float (/ ix x) 1.0d0)))
1105        (setq ra (+ (%rational-log x 1.0d0)
1106                    (/ (log1+ (* r r)) 2)))
1107        (setq ia (atan r (if (minusp rx) -1.0d0 1.0d0)))))
1108    (if (typep prototype 'short-float)
1109      (complex (%short-float ra) (%short-float ia))
1110      (complex ra ia))))
1111
1112(defun exp (x)
1113  "Return e raised to the power NUMBER."
1114  (typecase x
1115    (complex (* (exp (realpart x)) (cis (imagpart x))))
1116    (double-float (%double-float-exp! x (%make-dfloat)))
1117    (t
1118     #+32-bit-target
1119     (target::with-stack-short-floats ((sx x))
1120       (%single-float-exp! sx (%make-sfloat)))
1121     #+64-bit-target
1122     (%single-float-exp (%short-float x)))))
1123
1124
1125(defun positive-realpart-p (n)
1126  (> (realpart n) 0))
1127
1128;;; (It may be possible to do something with rational exponents, e.g. so that
1129;;;       (expt x 1/2) == (sqrt x)
1130;;;       (expt x 3/2) == (expt (sqrt x) 3)      ** NOT (sqrt (expt x 3)) !! **
1131;;;                      or (* x (sqrt x))
1132;;;       (expt x 1/8) == (sqrt (sqrt (sqrt x)))
1133;;;    even, in the case of rational x, returning a rational result if possible.)
1134;;;
1135(defun expt (b e)
1136  "Return BASE raised to the POWER."
1137  (cond ((zerop e) (1+ (* b e)))
1138        ((integerp e)
1139         (if (minusp e) (/ 1 (%integer-power b (- e))) (%integer-power b e)))
1140        ((zerop b)
1141         (if (plusp (realpart e)) (* b e) (report-bad-arg e '(satisfies plusp))))
1142        ((and (realp b) (plusp b) (realp e)
1143              ; escape out very small or very large rationals
1144              ; - avoid problems converting to float
1145              (typecase b
1146                (bignum (<= b most-positive-short-float))
1147                (ratio (cond ((< b 0.5)
1148                              (>= b least-positive-normalized-short-float))
1149                             ((> b 3)
1150                              (<= b most-positive-short-float))))
1151                (t t)))
1152         ;; assumes that the library routines are accurate
1153         ;; (if not, just excise this whole clause)
1154         (if (or (typep b 'double-float)
1155                 (typep e 'double-float))
1156           (with-stack-double-floats ((b1 b)
1157                                      (e1 e))
1158             (%double-float-expt! b1 e1 (%make-dfloat)))
1159           #+32-bit-target
1160           (target::with-stack-short-floats ((b1 b)
1161                                             (e1 e))
1162             (%single-float-expt! b1 e1 (%make-sfloat)))
1163           #+64-bit-target
1164           (%single-float-expt (%short-float b) (%short-float e))
1165           ))
1166        ((typep b 'rational)
1167         (let ((r (exp (* e (%rational-log b 1.0d0)))))
1168           (cond ((typep (realpart e) 'double-float)
1169                  r)
1170                 ((typep r 'complex)
1171                  (complex (%short-float (realpart r)) (%short-float (imagpart r))))
1172                 (t
1173                  (%short-float r)))))
1174        ((typep (realpart b) 'rational)
1175         (let ((r (exp (* e (%rational-complex-log b 1.0d0)))))
1176           (if (typep (realpart e) 'double-float)
1177             r
1178             (complex (%short-float (realpart r)) (%short-float (imagpart r))))))
1179        (t
1180         ;; type upgrade b without losing -0.0 ...
1181         (let ((r (exp (* e (log-e (* b 1.0d0))))))
1182           (cond ((or (typep (realpart b) 'double-float)
1183                      (typep (realpart e) 'double-float))
1184                  r)
1185                 ((typep r 'complex)
1186                  (complex (%short-float (realpart r)) (%short-float (imagpart r))))
1187                 (t
1188                  (%short-float r)))))))
1189
1190
1191       
1192(defun sqrt (x &aux a b)
1193  "Return the square root of NUMBER."
1194  (cond ((zerop x) x)
1195        ((complexp x)
1196         (let ((rx (realpart x))
1197               (ix (imagpart x)))
1198           (cond ((rationalp rx)
1199                  (if (zerop rx)
1200                    (let ((s (sqrt (/ (abs ix) 2))))
1201                      (complex s (if (minusp ix) (- s) s)))
1202                    (let* ((s (+ (* rx rx) (* ix ix)))
1203                           (d (if (ratiop s)
1204                                (/ (isqrt (%numerator s))
1205                                   (isqrt (%denominator s)))
1206                                (isqrt s))))
1207                      (unless (eql s (* d d))
1208                        (setf d (%double-float-hypot (%double-float rx)
1209                                                     (%double-float ix))))
1210                      (cond ((minusp rx)
1211                             (setq b (sqrt (/ (- d rx) 2)))
1212                             (when (minusp ix)
1213                               (setq b (- b)))
1214                             (setq a (/ ix (+ b b))))
1215                            (t
1216                             (setq a (sqrt (/ (+ rx d) 2)))
1217                             (setq b (/ ix (+ a a)))))
1218                      (if (rationalp a)
1219                        (complex a b)
1220                        (complex (%short-float a) (%short-float b))))))
1221                 ((minusp rx)
1222                  (if (zerop ix)
1223                    (complex 0 (float-sign ix (sqrt (- rx))))
1224                    (let ((shift (cond ((< rx -1) -3)
1225                                       ((and (> rx -5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25)
1226                                       (t -1))))
1227                      (setq rx (scale-float rx shift))
1228                      (let ((s (fsqrt (- (abs (complex rx (scale-float ix shift))) rx))))
1229                        (setq b (scale-float s (ash (- -1 shift) -1)))
1230                        (when (minusp ix)
1231                          (setq b (- b)))
1232                        (setq a (/ ix (scale-float b 1)))
1233                        (complex a b)))))
1234                 (t
1235                  (if (zerop ix)
1236                    (complex (sqrt rx) ix)
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 (+ rx (abs (complex rx (scale-float ix shift)))))))
1242                        (setq a (scale-float s (ash (- -1 shift) -1)))
1243                        (setq b (/ ix (scale-float a 1)))
1244                        (complex a b))))))))
1245        ((minusp x) (complex 0 (sqrt (- x))))
1246        ((floatp x)
1247         (fsqrt x))
1248        ((and (integerp x) (eql x (* (setq a (isqrt x)) a))) a)
1249        ((and (ratiop x)
1250              (let ((n (numerator x))
1251                    d)
1252                (and (eql n (* (setq a (isqrt n)) a))
1253                     (eql (setq d (denominator x))
1254                          (* (setq b (isqrt d)) b)))))
1255         (/ a b))         
1256        (t
1257         (float (fsqrt (float x 0.0d0)) 1.0s0))))
1258
1259
1260
1261(defun asin (x)
1262  "Return the arc sine of NUMBER."
1263  (cond ((and (typep x 'double-float)
1264              (locally (declare (type double-float x))
1265                (and (<= -1.0d0 x)
1266                     (<= x 1.0d0))))
1267         (%double-float-asin! x (%make-dfloat)))
1268        ((and (typep x 'single-float)
1269              (locally (declare (type single-float x))
1270                (and (<= -1.0s0 x)
1271                     (<= x 1.0s0))))
1272         #+32-bit-target
1273         (%single-float-asin! x (%make-sfloat))
1274         #+64-bit-target
1275         (%single-float-asin x))
1276        ((and (typep x 'rational)
1277              (<= (abs x) 1))
1278         (if (integerp x)
1279           (case x
1280             (0 0.0s0)                          ; or simply 0 ??
1281             (1 single-float-half-pi)
1282             (-1 #.(- single-float-half-pi)))
1283           (atan x (sqrt (- 1 (* x x))))))
1284        (t
1285         (%complex-asin/acos x nil))
1286        ))
1287
1288
1289(defun acos (x)
1290  "Return the arc cosine of NUMBER."
1291  (cond ((and (typep x 'double-float)
1292              (locally (declare (type double-float x))
1293                (and (<= -1.0d0 x)
1294                     (<= x 1.0d0))))
1295         (%double-float-acos! x (%make-dfloat)))
1296        ((and (typep x 'short-float)
1297              (locally (declare (type short-float x))
1298                (and (<= -1.0s0 x)
1299                     (<= x 1.0s0))))
1300         #+32-bit-target
1301         (%single-float-acos! x (%make-sfloat))
1302         #+64-bit-target
1303         (%single-float-acos x))
1304        ((and (typep x 'rational)
1305              (<= (abs x) 1))
1306         (if (integerp x)
1307           (case x
1308             (0 single-float-half-pi)
1309             (1 0.0s0)                          ; or simply 0 ??
1310             (-1 single-float-pi))
1311           (atan (sqrt (- 1 (* x x))) x)))
1312        (t
1313         (%complex-asin/acos x t))
1314        ))
1315
1316;;; combined complex asin/acos routine
1317;;; argument acos is true for acos(z); false for asin(z)
1318;;;
1319;;; based on Hull, Fairgrieve & Tang, ACM TMS 23, 3, 299-335 (Sept. 1997)
1320(defun %complex-asin/acos (z acos)
1321  (let* ((rx (realpart z))
1322         (ix (imagpart z))
1323         (x (abs rx))
1324         (y (abs ix))
1325         (m (max x y)))
1326    (if (> m 1.8014399E+16)
1327      ;; Large argument escape
1328      (let ((log-s 0))
1329        (if (typep m 'double-float)
1330          (if (> m #.(/ most-positive-double-float 2))
1331            (setq log-s double-float-log2)
1332            (setq z (* 2 z)))
1333          (if (> m #.(/ most-positive-short-float 2))
1334            (setq log-s single-float-log2)
1335            (setq z (* 2 z))))
1336        (if acos
1337          (i* (+ log-s (log-e z))
1338              (if (minusp ix) +1 -1))
1339          (if (minusp ix)
1340            (i* (+ log-s (log-e (i* z 1))) -1)
1341            (i* (+ log-s (log-e (i* z -1))) 1))))
1342      (let ((qrx rx)
1343            (qx x)
1344            x-1 y2 s)
1345        (cond ((rationalp rx)
1346               (setq x-1 (float (abs (- x 1))))
1347               (setq rx (float rx))
1348               (setq x (abs rx))
1349               (setq y (float y))
1350               (setq y2 (* y y))
1351               (setq s (cond ((zerop x-1)
1352                              y)
1353                             ((> y x-1)
1354                              (let ((c (/ x-1 y)))
1355                                (* y (sqrt (1+ (* c c))))))
1356                             (t
1357                              (let ((c (/ y x-1)))
1358                                (* x-1 (sqrt (1+ (* c c)))))))))
1359              (t
1360               (setq x-1 (abs (- x 1)))
1361               (setq y2 (* y y))
1362               (setq s (if (zerop x-1)
1363                         y
1364                         (sqrt (+ (* x-1 x-1) y2))))))
1365        (let* ((x+1 (+ x 1))
1366               (r (sqrt (+ (* x+1 x+1) y2)))
1367               (a (/ (+ r s) 2))
1368               (b (/ rx a))
1369               (ra (if (<= (abs b) 0.6417)
1370                     (if acos (acos b) (asin b))
1371                     (let* ((r+x+1 (+ r x+1))
1372                            (s+x-1 (+ s x-1))
1373                            (a+x (+ a x))
1374                            (ry (if (<= qx 1)
1375                                  (let ((aa (+ (/ y2 r+x+1) s+x-1)))
1376                                    (sqrt (/ (* a+x aa) 2)))
1377                                  (let ((aa (+ (/ a+x r+x+1) (/ a+x s+x-1))))
1378                                    (* y (sqrt (/ aa 2)))))))
1379                       (if acos (atan ry rx) (atan rx ry)))))
1380               (ia (if (<= a 1.5)
1381                     (let* ((r+x+1 (+ r x+1))
1382                            (s+x-1 (+ s x-1))
1383                            (ll (if (< qx 1)
1384                                  (let* ((aa (/ (+ (/ 1 r+x+1) (/ 1 s+x-1)) 2)))
1385                                    (+ (* aa y2) (* y (sqrt (* aa (1+ a))))))
1386                                  (let* ((a-1 (/ (+ (/ y2 r+x+1) s+x-1) 2)))
1387                                    (+ a-1 (sqrt (* a-1 (1+ a))))))))
1388                       (log1+ ll))
1389                     (log (+ a (sqrt (1- (* a a))))))))
1390          ;; final fixup of signs
1391          (if acos
1392            (if (complexp z)
1393              (if (typep ix 'float)
1394                (setq ia (float-sign (- ix) ia))
1395                (if (plusp ix)
1396                  (setq ia (- ia))))
1397              (if (< qrx -1)
1398                (setq ia (- ia))))
1399            (if (complexp z)
1400              (if (typep ix 'float)
1401                (setq ia (float-sign ix ia))
1402                (if (minusp ix)
1403                  (setq ia (- ia))))
1404              (if (> qrx 1)
1405                (setq ia (- ia)))))
1406          (complex ra ia))))))
1407
1408
1409(defun fsqrt (x)
1410  (etypecase x
1411    (double-float (%double-float-sqrt! x (%make-dfloat)))
1412    (single-float
1413     #+32-bit-target
1414     (%single-float-sqrt! x (%make-sfloat))
1415     #+64-bit-target
1416     (%single-float-sqrt x))))
1417
1418
1419
1420(defun %df-atan2 (y x)
1421  (if (zerop x)
1422    (if (zerop y)
1423      (if (plusp (float-sign x))
1424        (if (eql y -0.0d0)
1425          -0.0d0
1426          0.0d0)
1427        (float-sign y pi))
1428      (float-sign y double-float-half-pi))
1429    (%double-float-atan2! y x (%make-dfloat))))
1430
1431#+32-bit-target
1432(defun %sf-atan2! (y x)
1433  (if (zerop x)
1434    (if (zerop y)
1435      (if (plusp (float-sign x))
1436        ;; Don't return Y (which may be stack-consed) here.
1437        ;; We know that (ZEROP Y) is true, so:
1438        (if (eql y -0.0s0)
1439          -0.0s0
1440          0.0s0)
1441        (float-sign y single-float-pi))
1442      (float-sign y single-float-half-pi))
1443    (%single-float-atan2! y x (%make-sfloat))))
1444
1445#+64-bit-target
1446(defun %sf-atan2 (y x)
1447  (if (zerop x)
1448    (if (zerop y)
1449      (if (plusp (float-sign x))
1450        y
1451        (float-sign y single-float-pi))
1452      (float-sign y single-float-half-pi))
1453    (%single-float-atan2 y x)))
1454
1455#+64-bit-target
1456(defun %short-float-exp (n)
1457  (let* ((bits (single-float-bits n)))
1458    (declare (type (unsigned-byte 32) bits))
1459    (ldb (byte IEEE-single-float-exponent-width IEEE-single-float-exponent-offset) bits)))
1460
1461
1462#+64-bit-target
1463(defun set-%short-float-exp (float exp)
1464  (host-single-float-from-unsigned-byte-32
1465   (dpb exp
1466        (byte IEEE-single-float-exponent-width
1467              IEEE-single-float-exponent-offset)
1468        (the (unsigned-byte 32) (single-float-bits float)))))
1469
1470#+64-bit-target
1471(defun %%scale-sfloat (float int)
1472  (* (the single-float float)
1473     (the single-float (host-single-float-from-unsigned-byte-32
1474                        (dpb int
1475                             (byte IEEE-single-float-exponent-width
1476                                   IEEE-single-float-exponent-offset)
1477                             0)))))
1478
1479#+64-bit-target
1480(defun %double-float-exp (n)
1481  (let* ((highword (double-float-bits n)))
1482    (declare (fixnum highword))
1483    (logand (1- (ash 1 IEEE-double-float-exponent-width))
1484            (ash highword (- (- IEEE-double-float-exponent-offset 32))))))
1485
1486#+64-bit-target
1487(defun set-%double-float-exp (float exp)
1488  (let* ((highword (double-float-bits float)))
1489    (declare (fixnum highword))
1490    (setf (uvref float target::double-float.val-high-cell)
1491          (dpb exp
1492               (byte IEEE-double-float-exponent-width
1493                     (- IEEE-double-float-exponent-offset 32))
1494               highword))
1495    exp))
1496
1497#+64-bit-target
1498(defun %integer-decode-double-float (f)
1499  (multiple-value-bind (hiword loword) (double-float-bits f)
1500    (declare (type (unsigned-byte 32) hiword loword))
1501    (let* ((exp (ldb (byte IEEE-double-float-exponent-width
1502                           (- IEEE-double-float-exponent-offset 32))
1503                     hiword))
1504           (mantissa (logior
1505                      (the fixnum
1506                        (dpb (ldb (byte (- IEEE-double-float-mantissa-width 32)
1507                                        IEEE-double-float-mantissa-offset)
1508                                  hiword)
1509                             (byte (- IEEE-double-float-mantissa-width 32)
1510                                   32)
1511                             loword))
1512                      (if (zerop exp)
1513                        0
1514                        (ash 1 IEEE-double-float-hidden-bit))))
1515           (sign (if (logbitp 31 hiword) -1 1)))
1516      (declare (fixnum exp mantissa sign))
1517      (values (ldb (byte 25 28) mantissa)
1518              (ldb (byte 28 0) mantissa)
1519              exp
1520              sign))))
1521
1522;;; end of l0-float.lisp
Note: See TracBrowser for help on using the repository browser.