source: trunk/abcl/src/org/armedbear/lisp/numbers.lisp @ 14207

Last change on this file since 14207 was 11391, checked in by vvoutilainen, 12 years ago

ABCL license is GPL + Classpath exception. This was intended
by Peter Graves, the original author. For reference, see
http://sourceforge.net/mailarchive/forum.php?thread_name=20040721115302.839%40prufrock&forum_name=armedbear-j-announce

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.0 KB
Line 
1;;; numbers.lisp
2;;;
3;;; Copyright (C) 2003-2006 Peter Graves
4;;; $Id: numbers.lisp 11391 2008-11-15 22:38:34Z vvoutilainen $
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;;; As a special exception, the copyright holders of this library give you
21;;; permission to link this library with independent modules to produce an
22;;; executable, regardless of the license terms of these independent
23;;; modules, and to copy and distribute the resulting executable under
24;;; terms of your choice, provided that you also meet, for each linked
25;;; independent module, the terms and conditions of the license of that
26;;; module.  An independent module is a module which is not derived from
27;;; or based on this library.  If you modify this library, you may extend
28;;; this exception to your version of the library, but you are not
29;;; obligated to do so.  If you do not wish to do so, delete this
30;;; exception statement from your version.
31
32;;; Adapted from CMUCL/SBCL.
33
34(in-package "SYSTEM")
35
36(defun signum (number)
37  "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
38  (if (zerop number)
39      number
40      (if (rationalp number)
41    (if (plusp number) 1 -1)
42    (/ number (abs number)))))
43
44(defun round (number &optional (divisor 1))
45  "Rounds number (or number/divisor) to nearest integer.
46   The second returned value is the remainder."
47  (multiple-value-bind (tru rem) (truncate number divisor)
48    (if (zerop rem)
49        (values tru rem)
50        (let ((thresh (/ (abs divisor) 2)))
51          (cond ((or (> rem thresh)
52                     (and (= rem thresh) (oddp tru)))
53                 (if (minusp divisor)
54                     (values (- tru 1) (+ rem divisor))
55                     (values (+ tru 1) (- rem divisor))))
56                ((let ((-thresh (- thresh)))
57                   (or (< rem -thresh)
58                       (and (= rem -thresh) (oddp tru))))
59                 (if (minusp divisor)
60                     (values (+ tru 1) (- rem divisor))
61                     (values (- tru 1) (+ rem divisor))))
62                (t (values tru rem)))))))
63
64(defun ffloor (number &optional (divisor 1))
65  "Same as FLOOR, but returns first value as a float."
66  (multiple-value-bind (tru rem) (ftruncate number divisor)
67    (if (and (not (zerop rem))
68       (if (minusp divisor)
69     (plusp number)
70     (minusp number)))
71  (values (1- tru) (+ rem divisor))
72  (values tru rem))))
73
74(defun fceiling (number &optional (divisor 1))
75  "Same as CEILING, but returns first value as a float."
76  (multiple-value-bind (tru rem) (ftruncate number divisor)
77    (if (and (not (zerop rem))
78       (if (minusp divisor)
79     (minusp number)
80     (plusp number)))
81  (values (+ tru 1) (- rem divisor))
82  (values tru rem))))
83
84(defun fround (number &optional (divisor 1))
85  "Same as ROUND, but returns first value as a float."
86  (multiple-value-bind (res rem)
87    (round number divisor)
88    (values (float res (if (floatp rem) rem 1.0)) rem)))
89
90;;; FIXME
91(defun rationalize (number)
92  (rational number))
93
94(defun gcd (&rest integers)
95  (cond ((null integers)
96         0)
97  ((null (cdr integers))
98         (let ((n (car integers)))
99           (if (integerp n)
100               (abs n)
101               (error 'type-error :datum n :expected-type 'integer))))
102  (t
103   (do ((gcd (car integers) (gcd-2 gcd (car rest)))
104        (rest (cdr integers) (cdr rest)))
105       ((null rest) gcd)))))
106
107;;; From discussion on comp.lang.lisp and Akira Kurihara.
108(defun isqrt (natural)
109  "Returns the root of the nearest integer less than natural which is a perfect
110   square."
111  (unless (and (integerp natural) (not (minusp natural)))
112    (error 'simple-type-error
113           :format-control "The value ~A is not a non-negative real number."
114           :format-arguments (list natural)))
115  (if (and (fixnump natural) (<= natural 24))
116      (cond ((> natural 15) 4)
117      ((> natural  8) 3)
118      ((> natural  3) 2)
119      ((> natural  0) 1)
120      (t 0))
121      (let* ((n-len-quarter (ash (integer-length natural) -2))
122       (n-half (ash natural (- (ash n-len-quarter 1))))
123       (n-half-isqrt (isqrt n-half))
124       (init-value (ash (1+ n-half-isqrt) n-len-quarter)))
125  (loop
126    (let ((iterated-value
127     (ash (+ init-value (truncate natural init-value)) -1)))
128      (unless (< iterated-value init-value)
129        (return init-value))
130      (setq init-value iterated-value))))))
131
132;; FIXME Need to add support for denormalized floats!
133
134;; "FLOAT-PRECISION returns the number of significant radix b digits present in
135;; FLOAT; if FLOAT is a float zero, then the result is an integer zero."
136
137;; "For normalized floats, the results of FLOAT-DIGITS and FLOAT-PRECISION are
138;; the same, but the precision is less than the number of representation digits
139;; for a denormalized or zero number.
140(defun float-precision (float)
141  (if (floatp float)
142      (cond ((zerop float)
143             0)
144            ((typep float 'single-float)
145             24)
146            ((typep float 'double-float)
147             53)
148            (t
149             ;; Shouldn't get here!
150             (aver nil)))
151      (error 'simple-type-error
152             :format-control "~S is not of type FLOAT."
153             :format-arguments (list float))))
154
155(defun decode-float (float)
156  (multiple-value-bind (significand exponent sign)
157      (integer-decode-float float)
158    (values (coerce (/ significand (expt 2 53)) 'float)
159            (+ exponent 53)
160            (if (minusp sign) -1.0 1.0))))
161
162(defun conjugate (number)
163  (etypecase number
164    (complex
165     (complex (realpart number) (- (imagpart number))))
166    (number
167     number)))
168
169(defun phase (number)
170  "Returns the angle part of the polar representation of a complex number.
171   For complex numbers, this is (atan (imagpart number) (realpart number)).
172   For non-complex positive numbers, this is 0.  For non-complex negative
173   numbers this is PI."
174  (etypecase number
175             (rational
176              (if (minusp number)
177                  (coerce pi 'single-float)
178                  0.0f0))
179             (single-float
180              (if (minusp (float-sign number))
181                  (coerce pi 'single-float)
182                  0.0f0))
183             (double-float
184              (if (minusp (float-sign number))
185                  (coerce pi 'double-float)
186                  0.0d0))
187             (complex
188              (if (zerop (realpart number))
189                  (coerce (* (/ pi 2) (signum (imagpart number)))
190                          (if (typep (imagpart number) 'double-float)
191                              'double-float 'single-float))
192                  (atan (imagpart number) (realpart number))))))
Note: See TracBrowser for help on using the repository browser.