source: trunk/j/src/org/armedbear/lisp/numbers.lisp @ 5728

Last change on this file since 5728 was 5728, checked in by piso, 17 years ago

jvmcompile => jvm-compile

File size: 8.7 KB
Line 
1;;; numbers.lisp
2;;;
3;;; Copyright (C) 2003-2004 Peter Graves
4;;; $Id: numbers.lisp,v 1.20 2004-02-09 13:07:21 piso Exp $
5;;;
6;;; This program is free software; you can redistribute it and/or
7;;; modify it under the terms of the GNU General Public License
8;;; as published by the Free Software Foundation; either version 2
9;;; of the License, or (at your option) any later version.
10;;;
11;;; This program is distributed in the hope that it will be useful,
12;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
13;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14;;; GNU General Public License for more details.
15;;;
16;;; You should have received a copy of the GNU General Public License
17;;; along with this program; if not, write to the Free Software
18;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
19
20;;; From CMUCL.
21
22(in-package "SYSTEM")
23
24(defun signum (number)
25  "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
26  (if (zerop number)
27      number
28      (if (rationalp number)
29    (if (plusp number) 1 -1)
30    (/ number (abs number)))))
31
32
33;;; If the numbers do not divide exactly and the result of (/ number divisor)
34;;; would be negative then decrement the quotient and augment the remainder by
35;;; the divisor.
36;;;
37(defun floor (number &optional (divisor 1))
38  "Returns the greatest integer not greater than number, or number/divisor.
39  The second returned value is (mod number divisor)."
40  (multiple-value-bind (tru rem) (truncate number divisor)
41    (if (and (not (zerop rem))
42       (if (minusp divisor)
43     (plusp number)
44     (minusp number)))
45  (values (1- tru) (+ rem divisor))
46  (values tru rem))))
47
48
49;;; If the numbers do not divide exactly and the result of (/ number divisor)
50;;; would be positive then increment the quotient and decrement the remainder by
51;;; the divisor.
52;;;
53(defun ceiling (number &optional (divisor 1))
54  "Returns the smallest integer not less than number, or number/divisor.
55  The second returned value is the remainder."
56  (multiple-value-bind (tru rem) (truncate number divisor)
57    (if (and (not (zerop rem))
58       (if (minusp divisor)
59     (minusp number)
60     (plusp number)))
61  (values (+ tru 1) (- rem divisor))
62  (values tru rem))))
63
64
65(defun round (number &optional (divisor 1))
66  "Rounds number (or number/divisor) to nearest integer.
67   The second returned value is the remainder."
68  (multiple-value-bind (tru rem) (truncate number divisor)
69    (if (zerop rem)
70        (values tru rem)
71        (let ((thresh (/ (abs divisor) 2)))
72          (cond ((or (> rem thresh)
73                     (and (= rem thresh) (oddp tru)))
74                 (if (minusp divisor)
75                     (values (- tru 1) (+ rem divisor))
76                     (values (+ tru 1) (- rem divisor))))
77                ((let ((-thresh (- thresh)))
78                   (or (< rem -thresh)
79                       (and (= rem -thresh) (oddp tru))))
80                 (if (minusp divisor)
81                     (values (+ tru 1) (- rem divisor))
82                     (values (- tru 1) (+ rem divisor))))
83                (t (values tru rem)))))))
84
85
86(defun rem (number divisor)
87  "Returns second result of TRUNCATE."
88  (multiple-value-bind (tru rem) (truncate number divisor)
89    (declare (ignore tru))
90    rem))
91
92
93(defun mod (number divisor)
94  "Returns second result of FLOOR."
95  (let ((rem (rem number divisor)))
96    (if (and (not (zerop rem))
97       (if (minusp divisor)
98     (plusp number)
99     (minusp number)))
100  (+ rem divisor)
101  rem)))
102
103
104(defun ftruncate (number &optional (divisor 1))
105  (multiple-value-bind (tru rem) (truncate number divisor)
106    (values (float tru) rem)))
107
108
109(defun ffloor (number &optional (divisor 1))
110  "Same as FLOOR, but returns first value as a float."
111  (multiple-value-bind (tru rem) (ftruncate number divisor)
112    (if (and (not (zerop rem))
113       (if (minusp divisor)
114     (plusp number)
115     (minusp number)))
116  (values (1- tru) (+ rem divisor))
117  (values tru rem))))
118
119
120(defun fceiling (number &optional (divisor 1))
121  "Same as CEILING, but returns first value as a float."
122  (multiple-value-bind (tru rem) (ftruncate number divisor)
123    (if (and (not (zerop rem))
124       (if (minusp divisor)
125     (minusp number)
126     (plusp number)))
127  (values (+ tru 1) (- rem divisor))
128  (values tru rem))))
129
130
131(defun fround (number &optional (divisor 1))
132  "Same as ROUND, but returns first value as a float."
133  (multiple-value-bind (res rem)
134    (round number divisor)
135    (values (float res (if (floatp rem) rem 1.0)) rem)))
136
137
138(defun rational (x)
139  "RATIONAL produces a rational number for any real numeric argument.  This is
140   more efficient than RATIONALIZE, but it assumes that floating-point is
141   completely accurate, giving a result that isn't as pretty."
142  (cond ((floatp x)
143         (multiple-value-bind (bits exp)
144           (integer-decode-float x)
145           (if (eql bits 0)
146               0
147               (let* ((int (if (minusp x) (- bits) bits))
148                      (digits (float-digits x))
149                      (ex (+ exp digits)))
150                 (if (minusp ex)
151                     (/ int (ash 1 (+ digits (- ex))))
152                     (/ (ash int ex) (ash 1 digits)))))))
153        ((rationalp x)
154         x)
155        (t
156         (error 'type-error "wrong type: ~S is not a real number" x))))
157
158
159;;; FIXME
160(defun rationalize (x)
161  (rational x))
162
163(defun gcd (&rest numbers)
164  "Returns the greatest common divisor of the arguments, which must be
165   integers.  Gcd with no arguments is defined to be 0."
166  (unless (every #'integerp numbers)
167    (error 'type-error))
168  (cond ((null numbers) 0)
169  ((null (cdr numbers)) (abs (car numbers)))
170  (t
171   (do ((gcd (car numbers)
172       (gcd-2 gcd (car rest)))
173        (rest (cdr numbers) (cdr rest)))
174       ((null rest) gcd)))))
175
176;;; From discussion on comp.lang.lisp and Akira Kurihara.
177(defun isqrt (n)
178  "Returns the root of the nearest integer less than n which is a perfect
179   square."
180  (declare (type unsigned-byte n) (values unsigned-byte))
181  (unless (and (numberp n) (not (minusp n)))
182    (error 'type-error "wrong type: ~A is not a non-negative real number" n))
183  ;; theoretically (> n 7) ,i.e., n-len-quarter > 0
184  (if (and (fixnump n) (<= n 24))
185      (cond ((> n 15) 4)
186      ((> n  8) 3)
187      ((> n  3) 2)
188      ((> n  0) 1)
189      (t 0))
190      (let* ((n-len-quarter (ash (integer-length n) -2))
191       (n-half (ash n (- (ash n-len-quarter 1))))
192       (n-half-isqrt (isqrt n-half))
193       (init-value (ash (1+ n-half-isqrt) n-len-quarter)))
194  (loop
195    (let ((iterated-value
196     (ash (+ init-value (truncate n init-value)) -1)))
197      (unless (< iterated-value init-value)
198        (return init-value))
199      (setq init-value iterated-value))))))
200
201
202(defun float-sign (float1 &optional (float2 (float 1 float1)))
203  "Returns a floating-point number that has the same sign as
204   float1 and, if float2 is given, has the same absolute value
205   as float2."
206  (* (if (minusp float1)
207   (float -1 float1)
208   (float 1 float1))
209     (abs float2)))
210
211(defun decode-float (float)
212  (multiple-value-bind (significand exponent sign) (integer-decode-float float)
213    (values (coerce (/ significand (expt 2 53)) 'float)
214            (+ exponent 53)
215            (if (minusp sign) -1.0 1.0))))
216
217(defun conjugate (number)
218  (etypecase number
219    (complex
220     (complex (realpart number) (- (imagpart number))))
221    (number
222     number)))
223
224(defun phase (number)
225  "Returns the angle part of the polar representation of a complex number.
226   For complex numbers, this is (atan (imagpart number) (realpart number)).
227   For non-complex positive numbers, this is 0.  For non-complex negative
228   numbers this is PI."
229  (etypecase number
230             (rational
231              (if (minusp number)
232                  (coerce pi 'single-float)
233                  0.0f0))
234             (single-float
235              (if (minusp (float-sign number))
236                  (coerce pi 'single-float)
237                  0.0f0))
238             (double-float
239              (if (minusp (float-sign number))
240                  (coerce pi 'double-float)
241                  0.0d0))
242             (complex
243              (if (zerop (realpart number))
244                  (* (/ pi 2) (signum (imagpart number)))
245                  (atan (imagpart number) (realpart number))))))
246
247(defun cis (theta)
248  "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
249  (if (complexp theta)
250      (error "argument to CIS is complex: ~S" theta)
251      (complex (cos theta) (sin theta))))
252
253(when (and (fboundp 'jvm::jvm-compile) (not (autoloadp 'jvm::jvm-compile)))
254  (mapcar #'jvm::jvm-compile '(floor
255                               ceiling
256                               round
257                               rem
258                               ftruncate
259                               ffloor
260                               fceiling
261                               fround
262                               rational
263                               rationalize)))
Note: See TracBrowser for help on using the repository browser.