← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/lsp/numlib.lsp

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2006-05-17 02:46:26 UTC
  • Revision ID: james.westby@ubuntu.com-20060517024626-lljr08ftv9g9vefl
Tags: upstream-0.9h-20060510
Import upstream version 0.9h-20060510

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lsp/numlib.lsp
 
1
;;;;  Copyright (c) 1984, Taiichi Yuasa and Masami Hagiya.
 
2
;;;;  Copyright (c) 1990, Giuseppe Attardi.
 
3
;;;;
 
4
;;;;    This program is free software; you can redistribute it and/or
 
5
;;;;    modify it under the terms of the GNU Library General Public
 
6
;;;;    License as published by the Free Software Foundation; either
 
7
;;;;    version 2 of the License, or (at your option) any later version.
 
8
;;;;
 
9
;;;;    See file '../Copyright' for full details.
 
10
 
 
11
;;;;                           number routines
 
12
 
 
13
(in-package "SYSTEM")
 
14
 
 
15
#-ecl-min
 
16
(ffi:clines "#include <math.h>")
 
17
 
 
18
#.
 
19
(flet ((binary-search (f min max)
 
20
         (do ((new (/ (+ min max) 2) (/ (+ min max) 2)))
 
21
             ((>= min max)
 
22
              max)
 
23
           (if (funcall f new)
 
24
               (if (= new max)
 
25
                   (return max)
 
26
                   (setq max new))
 
27
               (if (= new min)
 
28
                   (return max)
 
29
                   (setq min new)))))
 
30
       (epsilon+ (x)
 
31
         (/= (float 1 x) (+ (float 1 x) x)))
 
32
       (epsilon- (x)
 
33
         (/= (float 1 x) (- (float 1 x) x))))
 
34
  `(eval-when (compile load eval)
 
35
    (defconstant short-float-epsilon
 
36
      ,(binary-search #'epsilon+ (coerce 0 'short-float) (coerce 1 'short-float)))
 
37
    (defconstant single-float-epsilon
 
38
      ,(binary-search #'epsilon+ (coerce 0 'single-float) (coerce 1 'single-float)))
 
39
    (defconstant long-float-epsilon
 
40
      ,(binary-search #'epsilon+ (coerce 0 'double-float) (coerce 1 'double-float)))
 
41
    (defconstant double-float-epsilon
 
42
      ,(binary-search #'epsilon+ (coerce 0 'long-float) (coerce 1 'long-float)))
 
43
    (defconstant short-float-negative-epsilon
 
44
      ,(binary-search #'epsilon- (coerce 0 'short-float) (coerce 1 'short-float)))
 
45
    (defconstant single-float-negative-epsilon
 
46
      ,(binary-search #'epsilon- (coerce 0 'single-float) (coerce 1 'single-float)))
 
47
    (defconstant long-float-negative-epsilon
 
48
      ,(binary-search #'epsilon- (coerce 0 'double-float) (coerce 1 'double-float)))
 
49
    (defconstant double-float-negative-epsilon
 
50
      ,(binary-search #'epsilon- (coerce 0 'long-float) (coerce 1 'long-float)))))
 
51
 
 
52
(defconstant imag-one #C(0.0 1.0))
 
53
 
 
54
(defun isqrt (i)
 
55
  "Args: (integer)
 
56
Returns the integer square root of INTEGER."
 
57
       (unless (and (integerp i) (>= i 0))
 
58
               (error 'type-error :datum i :expected-type 'unsigned-byte))
 
59
       (if (zerop i)
 
60
           0
 
61
           (let ((n (integer-length i)))
 
62
                (do ((x (ash 1 (ceiling n 2)))
 
63
                     (y))
 
64
                    (nil)
 
65
                    (setq y (floor i x))
 
66
                    (when (<= x y)
 
67
                          (return x))
 
68
                    (setq x (floor (+ x y) 2))))))
 
69
 
 
70
(defun abs (z)
 
71
  "Args: (number)
 
72
Returns the absolute value of NUMBER."
 
73
  (if (complexp z)
 
74
      ;; Compute sqrt(x*x + y*y) carefully to prevent overflow.
 
75
      ;; Assume |x| >= |y|. Then sqrt(x*x + y*y) = |x|*sqrt(1 +(y/x)^2).
 
76
      (let* ((x (abs (realpart z)))
 
77
             (y (abs (imagpart z))))
 
78
        ;; Swap x and y so that |x| >= |y|.
 
79
        (if (< x y)
 
80
            (rotatef x y))
 
81
        (if (zerop x)
 
82
            x
 
83
            (let ((r (/ y x)))
 
84
              (* x (sqrt (+ 1 (* r r)))))))
 
85
      (if (minusp z)
 
86
          (- z)
 
87
          z)))
 
88
 
 
89
(defun phase (x)
 
90
  "Args: (number)
 
91
Returns the angle part (in radians) of the polar representation of NUMBER.
 
92
Returns zero for non-complex numbers."
 
93
  (if (zerop x)
 
94
    (if (eq x 0) 0.0 (float 0 (realpart x)))
 
95
    (atan (imagpart x) (realpart x))))
 
96
 
 
97
(defun signum (x)
 
98
  "Args: (number)
 
99
Returns a number that represents the sign of NUMBER.  Returns NUMBER If it is
 
100
zero.  Otherwise, returns the value of (/ NUMBER (ABS NUMBER))"
 
101
  (if (zerop x) x (/ x (abs x))))
 
102
 
 
103
(defun cis (x)
 
104
  "Args: (radians)
 
105
Returns a complex number whose realpart and imagpart are the values of (COS
 
106
RADIANS) and (SIN RADIANS) respectively."
 
107
  (exp (* imag-one x)))
 
108
 
 
109
(defun asin (x)
 
110
  "Args: (number)
 
111
Returns the arc sine of NUMBER."
 
112
  (if #+ecl-min t #-ecl-min (complexp x)
 
113
      (complex-asin x)
 
114
      #-ecl-min
 
115
      (let* ((x (float x))
 
116
             (xr (float x 1d0)))
 
117
        (declare (double-float xr))
 
118
        (if (and (<= -1.0 xr) (<= xr 1.0))
 
119
            (float (ffi::c-inline (xr) (:double) :double "asin(#0)" :one-liner t)
 
120
                   x)
 
121
            (complex-asin x)))))
 
122
 
 
123
;; Ported from CMUCL
 
124
(defun complex-asin (z)
 
125
  (declare (number z)
 
126
           (si::c-local))
 
127
  (let ((sqrt-1-z (sqrt (- 1 z)))
 
128
        (sqrt-1+z (sqrt (+ 1 z))))
 
129
    (complex (atan (realpart z) (realpart (* sqrt-1-z sqrt-1+z)))
 
130
             (asinh (imagpart (* (conjugate sqrt-1-z)
 
131
                                 sqrt-1+z))))))
 
132
 
 
133
(defun acos (x)
 
134
  "Args: (number)
 
135
Returns the arc cosine of NUMBER."
 
136
  (if #+ecl-min t #-ecl-min (complexp x)
 
137
      (complex-acos x)
 
138
      #-ecl-min
 
139
      (let* ((x (float x))
 
140
             (xr (float x 1d0)))
 
141
        (declare (double-float xr))
 
142
        (if (and (<= -1.0 xr) (<= xr 1.0))
 
143
            (float (ffi::c-inline (xr) (:double) :double "acos(#0)" :one-liner t)
 
144
                   (float x))
 
145
            (complex-acos x)))))
 
146
 
 
147
;; Ported from CMUCL
 
148
(defun complex-acos (z)
 
149
  (declare (number z)
 
150
           (si::c-local))
 
151
  (let ((sqrt-1+z (sqrt (+ 1 z)))
 
152
        (sqrt-1-z (sqrt (- 1 z))))
 
153
    (complex (* 2 (atan (realpart sqrt-1-z) (realpart sqrt-1+z)))
 
154
             (asinh (imagpart (* (conjugate sqrt-1+z)
 
155
                                 sqrt-1-z))))))
 
156
 
 
157
#+(and (not ecl-min) win32)
 
158
(progn
 
159
  (ffi:clines "double asinh(double x) { return log(x+sqrt(1.0+x*x)); }")
 
160
  (ffi:clines "double acosh(double x) { return log(x+sqrt((x-1)*(x+1))); }")
 
161
  (ffi:clines "double atanh(double x) { return log((1+x)/(1-x))/2; }"))
 
162
 
 
163
;; Ported from CMUCL
 
164
(defun asinh (x)
 
165
  "Args: (number)
 
166
Returns the hyperbolic arc sine of NUMBER."
 
167
  ;(log (+ x (sqrt (+ 1.0 (* x x)))))
 
168
  (if #+(or ecl-min) t #-(or ecl-min) (complexp x)
 
169
      (let* ((iz (complex (- (imagpart x)) (realpart x)))
 
170
             (result (complex-asin iz)))
 
171
        (complex (imagpart result)
 
172
                 (- (realpart result))))
 
173
      #-(or ecl-min)
 
174
      (float (ffi:c-inline (x) (:double) :double "asinh(#0)" :one-liner t)
 
175
             (float x))))
 
176
 
 
177
;; Ported from CMUCL
 
178
(defun acosh (x)
 
179
  "Args: (number)
 
180
Returns the hyperbolic arc cosine of NUMBER."
 
181
  ;(log (+ x (sqrt (* (1- x) (1+ x)))))
 
182
  (if #+(or ecl-min) t #-(or ecl-min) (complexp x)
 
183
      (complex-acos x)
 
184
      #-(or ecl-min)
 
185
      (let* ((x (float x))
 
186
             (xr (float x 1d0)))
 
187
        (declare (double-float xr))
 
188
        (if (<= 1.0 xr)
 
189
            (float (ffi::c-inline (xr) (:double) :double "acosh(#0)" :one-liner t)
 
190
                   (float x))
 
191
            (complex-acosh x)))))
 
192
 
 
193
(defun complex-acosh (z)
 
194
  (declare (number z) (si::c-local))
 
195
  (let ((sqrt-z-1 (sqrt (- z 1)))
 
196
        (sqrt-z+1 (sqrt (+ z 1))))
 
197
    (complex (asinh (realpart (* (conjugate sqrt-z-1)
 
198
                                 sqrt-z+1)))
 
199
             (* 2 (atan (imagpart sqrt-z-1) (realpart sqrt-z+1))))))
 
200
 
 
201
(defun atanh (x)
 
202
  "Args: (number)
 
203
Returns the hyperbolic arc tangent of NUMBER."
 
204
  ;(/ (- (log (1+ x)) (log (- 1 x))) 2)
 
205
  (if #+(or ecl-min) t #-(or ecl-min) (complexp x)
 
206
      (complex-atanh x)
 
207
      #-(or ecl-min)
 
208
      (let* ((x (float x))
 
209
             (xr (float x 1d0)))
 
210
        (declare (double-float xr))
 
211
        (if (and (<= -1.0 xr) (<= xr 1.0))
 
212
            (float (ffi::c-inline (xr) (:double) :double "atanh(#0)" :one-liner t)
 
213
                   (float x))
 
214
            (complex-atanh x)))))
 
215
 
 
216
(defun complex-atanh (z)
 
217
  (declare (number x) (si::c-local))
 
218
  (/ (- (log (1+ z)) (log (- 1 z))) 2))
 
219
 
 
220
(defun ffloor (x &optional (y 1))
 
221
  "Args: (number &optional (divisor 1))
 
222
Same as FLOOR, but returns a float as the first value."
 
223
  (multiple-value-bind (i r) (floor x y)
 
224
    (values (if (floatp r) (float i r) (float i)) r)))
 
225
 
 
226
(defun fceiling (x &optional (y 1.0s0))
 
227
  "Args: (number &optional (divisor 1))
 
228
Same as CEILING, but returns a float as the first value."
 
229
  (multiple-value-bind (i r) (ceiling x y)
 
230
    (values (if (floatp r) (float i r) (float i)) r)))
 
231
 
 
232
(defun ftruncate (x &optional (y 1.0s0))
 
233
  "Args: (number &optional (divisor 1))
 
234
Same as TRUNCATE, but returns a float as the first value."
 
235
  (multiple-value-bind (i r) (truncate x y)
 
236
    (values (if (floatp r) (float i r) (float i)) r)))
 
237
 
 
238
(defun fround (x &optional (y 1.0s0))
 
239
  "Args: (number &optional (divisor 1))
 
240
Same as ROUND, but returns a float as the first value."
 
241
  (multiple-value-bind (i r) (round x y)
 
242
    (values (if (floatp r) (float i r) (float i)) r)))
 
243
 
 
244
(defun logtest (x y)
 
245
  "Args: (integer1 integer2)
 
246
Equivalent to (NOT (ZEROP (LOGAND INTEGER1 INTEGER2)))."
 
247
  (not (zerop (logand x y))))
 
248
 
 
249
 
 
250
(defun byte (size position)
 
251
  "Args: (size position)
 
252
Returns a byte specifier of integers.  The value specifies the SIZE-bits byte
 
253
starting the least-significant-bit but POSITION bits of integers.  In ECL, a
 
254
byte specifier is represented by a dotted pair (SIZE . POSITION)."
 
255
  (cons size position))
 
256
 
 
257
(defun byte-size (bytespec)
 
258
  "Args: (byte)
 
259
Returns the size part (in ECL, the car part) of the byte specifier BYTE."
 
260
  (car bytespec))
 
261
 
 
262
(defun byte-position (bytespec)
 
263
  "Args: (byte)
 
264
Returns the position part (in ECL, the cdr part) of the byte specifier BYTE."
 
265
  (cdr bytespec))
 
266
 
 
267
(defun ldb (bytespec integer)
 
268
  "Args: (bytespec integer)
 
269
Extracts a byte from INTEGER at the specified byte position, right-justifies
 
270
the byte, and returns the result as an integer."
 
271
  (logandc2 (ash integer (- (byte-position bytespec)))
 
272
            (- (ash 1 (byte-size bytespec)))))
 
273
 
 
274
(defun ldb-test (bytespec integer)
 
275
  "Args: (bytespec integer)
 
276
Returns T if at least one bit of the specified byte is 1; NIL otherwise."
 
277
  (not (zerop (ldb bytespec integer))))
 
278
 
 
279
(defun mask-field (bytespec integer)
 
280
  "Args: (bytespec integer)
 
281
Extracts the specified byte from INTEGER and returns the result as an integer."
 
282
  (ash (ldb bytespec integer) (byte-position bytespec)))
 
283
 
 
284
(defun dpb (newbyte bytespec integer)
 
285
  "Args: (newbyte bytespec integer)
 
286
Replaces the specified byte of INTEGER with NEWBYTE (an integer) and returns
 
287
the result."
 
288
  (logxor integer
 
289
          (mask-field bytespec integer)
 
290
          (ash (logandc2 newbyte
 
291
                         (- (ash 1 (byte-size bytespec))))
 
292
               (byte-position bytespec))))
 
293
 
 
294
(defun deposit-field (newbyte bytespec integer)
 
295
  "Args: (integer1 bytespec integer2)
 
296
Returns an integer represented by the bit sequence obtained by replacing the
 
297
specified bits of INTEGER2 with the specified bits of INTEGER1."
 
298
  (dpb (ash newbyte (- (byte-position bytespec))) bytespec integer))