← Back to branch summary

~ubuntu-branches/ubuntu/quantal/gclcvs/quantal

« back to all changes in this revision

Viewing changes to lsp/gcl_numlib.lsp

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-06-24 15:13:46 UTC
  • Revision ID: james.westby@ubuntu.com-20040624151346-xh0xaaktyyp7aorc
Tags: 2.7.0-26
C_GC_OFFSET is 2 on m68k-linux

Show diffs side-by-side

added added

removed removed

Lines of Context:
lsp/gcl_numlib.lsp
 
1
;; Copyright (C) 1994 M. Hagiya, W. Schelter, T. Yuasa
 
2
 
 
3
;; This file is part of GNU Common Lisp, herein referred to as GCL
 
4
;;
 
5
;; GCL is free software; you can redistribute it and/or modify it under
 
6
;;  the terms of the GNU LIBRARY GENERAL PUBLIC LICENSE as published by
 
7
;; the Free Software Foundation; either version 2, or (at your option)
 
8
;; any later version.
 
9
;; 
 
10
;; GCL is distributed in the hope that it will be useful, but WITHOUT
 
11
;; ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 
12
;; FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public 
 
13
;; License for more details.
 
14
;; 
 
15
;; You should have received a copy of the GNU Library General Public License 
 
16
;; along with GCL; see the file COPYING.  If not, write to the Free Software
 
17
;; Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
 
18
 
 
19
 
 
20
;;;;    numlib.lsp
 
21
;;;;
 
22
;;;;                           number routines
 
23
 
 
24
 
 
25
(in-package 'lisp)
 
26
(export
 
27
 '(isqrt abs phase signum cis asin acos sinh cosh tanh
 
28
   asinh acosh atanh
 
29
   rational rationalize
 
30
   ffloor fround ftruncate fceiling
 
31
   lognand lognor logandc1 logandc2 logorc1 logorc2
 
32
   lognot logtest
 
33
   byte byte-size byte-position
 
34
   ldb ldb-test mask-field dpb deposit-field
 
35
   ))
 
36
 
 
37
 
 
38
(in-package 'system)
 
39
 
 
40
 
 
41
(proclaim '(optimize (safety 2) (space 3)))
 
42
 
 
43
 
 
44
(defconstant imag-one #C(0.0d0 1.0d0))
 
45
 
 
46
 
 
47
(defun isqrt (i)
 
48
       (unless (and (integerp i) (>= i 0))
 
49
               (error "~S is not a non-negative integer." i))
 
50
       (if (zerop i)
 
51
           0
 
52
           (let ((n (integer-length i)))
 
53
                (do ((x (ash 1 (ceiling n 2)))
 
54
                     (y))
 
55
                    (nil)
 
56
                    (setq y (floor i x))
 
57
                    (when (<= x y)
 
58
                          (return x))
 
59
                    (setq x (floor (+ x y) 2))))))
 
60
 
 
61
 
 
62
(defun abs (z)
 
63
  (cond ((complexp z)
 
64
         ;; Compute (sqrt (+ (* x x) (* y y))) carefully to prevent
 
65
         ;; overflow!
 
66
         (let* ((x (abs (realpart z)))
 
67
                (y (abs (imagpart z))))
 
68
           (if (< x y)
 
69
               (rotatef x y))
 
70
           (if (zerop x)
 
71
               x
 
72
               (let ((r (/  y x)))
 
73
                 (* x (sqrt (+ 1 (* r r))))))))
 
74
        (t                              ; Should this be (realp z) instead of t?
 
75
         (if (minusp z)
 
76
             (- z)
 
77
             z))))
 
78
 
 
79
 
 
80
(defun phase (x)
 
81
       (atan (imagpart x) (realpart x)))
 
82
 
 
83
(defun signum (x) (if (zerop x) x (/ x (abs x))))
 
84
 
 
85
(defun cis (x) (exp (* imag-one x)))
 
86
 
 
87
(defun asin (x)
 
88
       (let ((c (- (* imag-one
 
89
                      (log (+ (* imag-one x)
 
90
                              (sqrt (- 1.0d0 (* x x)))))))))
 
91
            (if (or (and (not (complexp x))
 
92
                          (<= x 1.0d0)
 
93
                           (>= x -1.0d0)
 
94
                            )
 
95
                        (zerop (imagpart c)))
 
96
                (realpart c)
 
97
                c)))
 
98
 
 
99
(defun acos (x)
 
100
       (let ((c (- (* imag-one
 
101
                      (log (+ x (* imag-one
 
102
                                   (sqrt (- 1.0d0 (* x x))))))))))
 
103
            (if (or (and (not (complexp x))
 
104
                          (<= x 1.0d0)
 
105
                           (>= x -1.0d0)
 
106
                            )
 
107
                        (zerop (imagpart c)))
 
108
                (realpart c)
 
109
                c)))
 
110
 
 
111
 
 
112
(defun sinh (z)
 
113
  (cond ((complexp z)
 
114
         ;; For complex Z, compute the real and imaginary parts
 
115
         ;; separately to get better precision.
 
116
         (let ((x (realpart z))
 
117
               (y (imagpart z)))
 
118
           (complex (* (sinh x) (cos y))
 
119
                    (* (cosh x) (sin y)))))
 
120
        (t                              ; Should this be (realp z) instead of t?
 
121
         (let ((limit #.(expt (* double-float-epsilon 45/2) 1/5)))
 
122
           (if (< (- limit) z limit)
 
123
               ;; For this region, write use the fact that sinh z =
 
124
               ;; z*exp(z)*[(1 - exp(-2z))/(2z)].  Then use the first
 
125
               ;; 4 terms in the Taylor series expansion of
 
126
               ;; (1-exp(-2z))/2/z.  series expansion of (1 -
 
127
               ;; exp(2*x)).  This is needed because there is severe
 
128
               ;; roundoff error calculating (1 - exp(-2z)) for z near
 
129
               ;; 0.
 
130
               (* z (exp z)
 
131
                  (- 1 (* z
 
132
                          (- 1 (* z
 
133
                                  (- 2/3 (* z
 
134
                                            (- 1/3 (* 2/15 z)))))))))
 
135
               (let ((e (exp z)))
 
136
                 (* 1/2 (- e (/ e)))))))))
 
137
 
 
138
;(defun sinh (x) (/ (- (exp x) (exp (- x))) 2.0d0))
 
139
 
 
140
(defun cosh (z)
 
141
  (cond ((complexp z)
 
142
         ;; For complex Z, compute the real and imaginary parts
 
143
         ;; separately to get better precision.
 
144
         (let ((x (realpart z))
 
145
               (y (imagpart z)))
 
146
           (complex (* (cosh x) (cos y))
 
147
                    (* (sinh x) (sin y)))))
 
148
        (t                              ; Should this be (realp z) instead of t?
 
149
         ;; For real Z, there's no chance of round-off error, so
 
150
         ;; direct evaluation is ok.
 
151
         (let ((e (exp z)))
 
152
           (* 1/2 (+ e (/ e)))))))
 
153
;(defun cosh (x) (/ (+ (exp x) (exp (- x))) 2.0d0))
 
154
(defun tanh (x) (/ (sinh x) (cosh x)))
 
155
 
 
156
(defun asinh (x) (log (+ x (sqrt (+ 1.0d0 (* x x))))))
 
157
;(defun acosh (x)
 
158
;  (log (+ x
 
159
;         (* (1+ x)
 
160
;            (sqrt (/ (1- x) (1+ x)))))))
 
161
;(defun acosh (x)
 
162
;       (log (+ x
 
163
;              (sqrt (* (1- x) (1+ x))))))
 
164
(defun acosh (x)
 
165
  (* 2 (log (+ (sqrt (/ (1+ x) 2)) (sqrt (/ (1- x) 2))))))
 
166
(defun atanh (x)
 
167
       (when (or (= x 1.0d0) (= x -1.0d0))
 
168
             (error "The argument, ~s, is a logarithmic singularity.~
 
169
                    ~%Don't be foolish, GLS."
 
170
                    x))
 
171
       (log (/ (1+ x) (sqrt (- 1.0d0 (* x x))))))
 
172
 
 
173
 
 
174
(defun rational (x)
 
175
  (etypecase x
 
176
    (float        
 
177
      (multiple-value-bind (i e s) (integer-decode-float x)
 
178
                           (if (>= s 0)
 
179
                               (* i (expt (float-radix x) e))
 
180
                             (- (* i (expt (float-radix x) e))))))
 
181
    (rational x)))
 
182
 
 
183
 
 
184
(setf (symbol-function 'rationalize) (symbol-function 'rational))
 
185
 
 
186
;; although the following is correct code in that it approximates the
 
187
;; x to within eps, it does not preserve (eql (float (rationalize x) x) x)
 
188
;; since the test for eql is more strict than the float-epsilon
 
189
 
 
190
;;; Rationalize originally by Skef Wholey.
 
191
;;; Obtained from Daniel L. Weinreb.
 
192
;(defun rationalize (x)
 
193
;  (typecase x
 
194
;    (rational x)
 
195
;    (short-float (rationalize-float x short-float-epsilon 1.0s0))
 
196
;    (long-float (rationalize-float x long-float-epsilon 1.0d0))
 
197
;    (otherwise (error "~S is neither rational nor float." x))))
 
198
;
 
199
;(defun rationalize-float (x eps one)
 
200
;  (cond ((minusp x) (- (rationalize (- x))))
 
201
;        ((zerop x) 0)
 
202
;        (t (let ((y ())
 
203
;                 (a ()))
 
204
;             (do ((xx x (setq y (/ one
 
205
;                                   (- xx (float a x)))))
 
206
;                  (num (setq a (truncate x))
 
207
;                       (+ (* (setq a (truncate y)) num) onum))
 
208
;                  (den 1 (+ (* a den) oden))
 
209
;                  (onum 1 num)
 
210
;                  (oden 0 den))
 
211
;                 ((and (not (zerop den))
 
212
;                       (not (> (abs (/ (- x (/ (float num x)
 
213
;                                               (float den x)))
 
214
;                                       x))
 
215
;                               eps)))
 
216
;                  (/ num den)))))))
 
217
 
 
218
 
 
219
(defun ffloor (x &optional (y 1.0s0))
 
220
       (multiple-value-bind (i r) (floor (float x) (float y))
 
221
        (values (float i r) r)))
 
222
 
 
223
(defun fceiling (x &optional (y 1.0s0))
 
224
       (multiple-value-bind (i r) (ceiling (float x) (float y))
 
225
        (values (float i r) r)))
 
226
 
 
227
(defun ftruncate (x &optional (y 1.0s0))
 
228
       (multiple-value-bind (i r) (truncate (float x) (float y))
 
229
        (values (float i r) r)))
 
230
 
 
231
(defun fround (x &optional (y 1.0s0))
 
232
       (multiple-value-bind (i r) (round (float x) (float y))
 
233
        (values (float i r) r)))
 
234
 
 
235
 
 
236
(defun lognand (x y) (boole boole-nand x y))
 
237
(defun lognor (x y) (boole boole-nor x y))
 
238
(defun logandc1 (x y) (boole boole-andc1 x y))
 
239
(defun logandc2 (x y) (boole boole-andc2 x y))
 
240
(defun logorc1 (x y) (boole boole-orc1 x y))
 
241
(defun logorc2 (x y) (boole boole-orc2 x y))
 
242
 
 
243
(defun lognot (x) (logxor -1 x))
 
244
(defun logtest (x y) (not (zerop (logand x y))))
 
245
 
 
246
 
 
247
(defun byte (size position)
 
248
  (cons size position))
 
249
 
 
250
(defun byte-size (bytespec)
 
251
  (car bytespec))
 
252
 
 
253
(defun byte-position (bytespec)
 
254
  (cdr bytespec))
 
255
 
 
256
(defun ldb (bytespec integer)
 
257
  (logandc2 (ash integer (- (byte-position bytespec)))
 
258
            (- (ash 1 (byte-size bytespec)))))
 
259
 
 
260
(defun ldb-test (bytespec integer)
 
261
  (not (zerop (ldb bytespec integer))))
 
262
 
 
263
(defun mask-field (bytespec integer)
 
264
  (ash (ldb bytespec integer) (byte-position bytespec)))
 
265
 
 
266
(defun dpb (newbyte bytespec integer)
 
267
  (logxor integer
 
268
          (mask-field bytespec integer)
 
269
          (ash (logandc2 newbyte
 
270
                         (- (ash 1 (byte-size bytespec))))
 
271
               (byte-position bytespec))))
 
272
 
 
273
(defun deposit-field (newbyte bytespec integer)
 
274
  (dpb (ash newbyte (- (byte-position bytespec))) bytespec integer))