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

Last change on this file since 3557 was 3557, checked in by piso, 19 years ago

ISQRT

File size: 6.2 KB
Line 
1;;; numbers.lisp
2;;;
3;;; Copyright (C) 2003 Peter Graves
4;;; $Id: numbers.lisp,v 1.9 2003-09-02 17:49:54 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;;; If the numbers do not divide exactly and the result of (/ number divisor)
25;;; would be negative then decrement the quotient and augment the remainder by
26;;; the divisor.
27;;;
28(defun floor (number &optional (divisor 1))
29  "Returns the greatest integer not greater than number, or number/divisor.
30  The second returned value is (mod number divisor)."
31  (multiple-value-bind (tru rem) (truncate number divisor)
32    (if (and (not (zerop rem))
33       (if (minusp divisor)
34     (plusp number)
35     (minusp number)))
36  (values (1- tru) (+ rem divisor))
37  (values tru rem))))
38
39
40;;; If the numbers do not divide exactly and the result of (/ number divisor)
41;;; would be positive then increment the quotient and decrement the remainder by
42;;; the divisor.
43;;;
44(defun ceiling (number &optional (divisor 1))
45  "Returns the smallest integer not less than number, or number/divisor.
46  The second returned value is the remainder."
47  (multiple-value-bind (tru rem) (truncate number divisor)
48    (if (and (not (zerop rem))
49       (if (minusp divisor)
50     (minusp number)
51     (plusp number)))
52  (values (+ tru 1) (- rem divisor))
53  (values tru rem))))
54
55
56(defun round (number &optional (divisor 1))
57  "Rounds number (or number/divisor) to nearest integer.
58   The second returned value is the remainder."
59  (multiple-value-bind (tru rem) (truncate number divisor)
60    (if (zerop rem)
61        (values tru rem)
62        (let ((thresh (/ (abs divisor) 2)))
63          (cond ((or (> rem thresh)
64                     (and (= rem thresh) (oddp tru)))
65                 (if (minusp divisor)
66                     (values (- tru 1) (+ rem divisor))
67                     (values (+ tru 1) (- rem divisor))))
68                ((let ((-thresh (- thresh)))
69                   (or (< rem -thresh)
70                       (and (= rem -thresh) (oddp tru))))
71                 (if (minusp divisor)
72                     (values (+ tru 1) (- rem divisor))
73                     (values (- tru 1) (+ rem divisor))))
74                (t (values tru rem)))))))
75
76
77(defun rem (number divisor)
78  "Returns second result of TRUNCATE."
79  (multiple-value-bind (tru rem) (truncate number divisor)
80    (declare (ignore tru))
81    rem))
82
83
84(defun mod (number divisor)
85  "Returns second result of FLOOR."
86  (let ((rem (rem number divisor)))
87    (if (and (not (zerop rem))
88       (if (minusp divisor)
89     (plusp number)
90     (minusp number)))
91  (+ rem divisor)
92  rem)))
93
94
95(defun ftruncate (number &optional (divisor 1))
96  (multiple-value-bind (tru rem) (truncate number divisor)
97    (values (float tru) rem)))
98
99
100(defun ffloor (number &optional (divisor 1))
101  "Same as FLOOR, but returns first value as a float."
102  (multiple-value-bind (tru rem) (ftruncate number divisor)
103    (if (and (not (zerop rem))
104       (if (minusp divisor)
105     (plusp number)
106     (minusp number)))
107  (values (1- tru) (+ rem divisor))
108  (values tru rem))))
109
110
111(defun fceiling (number &optional (divisor 1))
112  "Same as CEILING, but returns first value as a float."
113  (multiple-value-bind (tru rem) (ftruncate number divisor)
114    (if (and (not (zerop rem))
115       (if (minusp divisor)
116     (minusp number)
117     (plusp number)))
118  (values (+ tru 1) (- rem divisor))
119  (values tru rem))))
120
121
122(defun fround (number &optional (divisor 1))
123  "Same as ROUND, but returns first value as a float."
124  (multiple-value-bind (res rem)
125    (round number divisor)
126    (values (float res (if (floatp rem) rem 1.0)) rem)))
127
128
129(defun rational (x)
130  "RATIONAL produces a rational number for any real numeric argument.  This is
131   more efficient than RATIONALIZE, but it assumes that floating-point is
132   completely accurate, giving a result that isn't as pretty."
133  (cond ((floatp x)
134         (multiple-value-bind (bits exp)
135           (integer-decode-float x)
136           (if (eql bits 0)
137               0
138               (let* ((int (if (minusp x) (- bits) bits))
139                      (digits (float-digits x))
140                      (ex (+ exp digits)))
141                 (if (minusp ex)
142                     (/ int (ash 1 (+ digits (- ex))))
143                     (/ (ash int ex) (ash 1 digits)))))))
144        ((rationalp x)
145         x)
146        (t
147         (error 'type-error "wrong type: ~S is not a real number" x))))
148
149;;; FIXME
150(defun rationalize (x)
151  (rational x))
152
153(when (and (find-package "JVM")
154           (fboundp 'jvm::jvm-compile))
155  (mapcar #'jvm::jvm-compile '(floor
156                               ceiling
157                               round
158                               rem
159                               ftruncate
160                               ffloor
161                               fceiling
162                               fround
163                               rational
164                               rationalize)))
165
166
167
168;;; From discussion on comp.lang.lisp and Akira Kurihara.
169(defun isqrt (n)
170  "Returns the root of the nearest integer less than n which is a perfect
171   square."
172  (declare (type unsigned-byte n) (values unsigned-byte))
173  ;; theoretically (> n 7) ,i.e., n-len-quarter > 0
174  (if (and (fixnump n) (<= n 24))
175      (cond ((> n 15) 4)
176      ((> n  8) 3)
177      ((> n  3) 2)
178      ((> n  0) 1)
179      (t 0))
180      (let* ((n-len-quarter (ash (integer-length n) -2))
181       (n-half (ash n (- (ash n-len-quarter 1))))
182       (n-half-isqrt (isqrt n-half))
183       (init-value (ash (1+ n-half-isqrt) n-len-quarter)))
184  (loop
185    (let ((iterated-value
186     (ash (+ init-value (truncate n init-value)) -1)))
187      (unless (< iterated-value init-value)
188        (return init-value))
189      (setq init-value iterated-value))))))
Note: See TracBrowser for help on using the repository browser.