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

Last change on this file since 8784 was 8784, checked in by piso, 16 years ago

SINGLE-FLOAT support.

File size: 5.8 KB
Line 
1;;; numbers.lisp
2;;;
3;;; Copyright (C) 2003-2005 Peter Graves
4;;; $Id: numbers.lisp,v 1.35 2005-03-17 15:04:23 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(defun round (number &optional (divisor 1))
33  "Rounds number (or number/divisor) to nearest integer.
34   The second returned value is the remainder."
35  (multiple-value-bind (tru rem) (truncate number divisor)
36    (if (zerop rem)
37        (values tru rem)
38        (let ((thresh (/ (abs divisor) 2)))
39          (cond ((or (> rem thresh)
40                     (and (= rem thresh) (oddp tru)))
41                 (if (minusp divisor)
42                     (values (- tru 1) (+ rem divisor))
43                     (values (+ tru 1) (- rem divisor))))
44                ((let ((-thresh (- thresh)))
45                   (or (< rem -thresh)
46                       (and (= rem -thresh) (oddp tru))))
47                 (if (minusp divisor)
48                     (values (+ tru 1) (- rem divisor))
49                     (values (- tru 1) (+ rem divisor))))
50                (t (values tru rem)))))))
51
52(defun ffloor (number &optional (divisor 1))
53  "Same as FLOOR, but returns first value as a float."
54  (multiple-value-bind (tru rem) (ftruncate number divisor)
55    (if (and (not (zerop rem))
56       (if (minusp divisor)
57     (plusp number)
58     (minusp number)))
59  (values (1- tru) (+ rem divisor))
60  (values tru rem))))
61
62(defun fceiling (number &optional (divisor 1))
63  "Same as CEILING, but returns first value as a float."
64  (multiple-value-bind (tru rem) (ftruncate number divisor)
65    (if (and (not (zerop rem))
66       (if (minusp divisor)
67     (minusp number)
68     (plusp number)))
69  (values (+ tru 1) (- rem divisor))
70  (values tru rem))))
71
72(defun fround (number &optional (divisor 1))
73  "Same as ROUND, but returns first value as a float."
74  (multiple-value-bind (res rem)
75    (round number divisor)
76    (values (float res (if (floatp rem) rem 1.0)) rem)))
77
78;;; FIXME
79(defun rationalize (number)
80  (rational number))
81
82(defun gcd (&rest integers)
83  "Returns the greatest common divisor of the arguments, which must be
84   integers.  Gcd with no arguments is defined to be 0."
85  (unless (every #'integerp integers)
86    (error 'type-error :datum (find-if-not #'integerp integers) :expected-type 'integer))
87  (cond ((null integers) 0)
88  ((null (cdr integers)) (abs (car integers)))
89  (t
90   (do ((gcd (car integers)
91       (gcd-2 gcd (car rest)))
92        (rest (cdr integers) (cdr rest)))
93       ((null rest) gcd)))))
94
95;;; From discussion on comp.lang.lisp and Akira Kurihara.
96(defun isqrt (natural)
97  "Returns the root of the nearest integer less than natural which is a perfect
98   square."
99  (unless (and (integerp natural) (not (minusp natural)))
100    (error 'simple-type-error
101           :format-control "The value ~A is not a non-negative real number."
102           :format-arguments (list natural)))
103  (if (and (fixnump natural) (<= natural 24))
104      (cond ((> natural 15) 4)
105      ((> natural  8) 3)
106      ((> natural  3) 2)
107      ((> natural  0) 1)
108      (t 0))
109      (let* ((n-len-quarter (ash (integer-length natural) -2))
110       (n-half (ash natural (- (ash n-len-quarter 1))))
111       (n-half-isqrt (isqrt n-half))
112       (init-value (ash (1+ n-half-isqrt) n-len-quarter)))
113  (loop
114    (let ((iterated-value
115     (ash (+ init-value (truncate natural init-value)) -1)))
116      (unless (< iterated-value init-value)
117        (return init-value))
118      (setq init-value iterated-value))))))
119
120(defun float-precision (float)
121  (if (floatp float)
122      (if (zerop float) 0 53)
123      (error 'simple-type-error
124             :format-control "~S is not of type FLOAT."
125             :format-arguments (list float))))
126
127(defun decode-float (float)
128  (multiple-value-bind (significand exponent sign) (integer-decode-float float)
129    (values (coerce (/ significand (expt 2 53)) 'float)
130            (+ exponent 53)
131            (if (minusp sign) -1.0 1.0))))
132
133(defun conjugate (number)
134  (etypecase number
135    (complex
136     (complex (realpart number) (- (imagpart number))))
137    (number
138     number)))
139
140(defun phase (number)
141  "Returns the angle part of the polar representation of a complex number.
142   For complex numbers, this is (atan (imagpart number) (realpart number)).
143   For non-complex positive numbers, this is 0.  For non-complex negative
144   numbers this is PI."
145  (etypecase number
146             (rational
147              (if (minusp number)
148                  (coerce pi 'single-float)
149                  0.0f0))
150             (single-float
151              (if (minusp (float-sign number))
152                  (coerce pi 'single-float)
153                  0.0f0))
154             (double-float
155              (if (minusp (float-sign number))
156                  (coerce pi 'double-float)
157                  0.0d0))
158             (complex
159              (if (zerop (realpart number))
160                  (coerce (* (/ pi 2) (signum (imagpart number)))
161                          (if (typep (imagpart number) 'double-float)
162                              'double-float 'single-float))
163                  (atan (imagpart number) (realpart number))))))
Note: See TracBrowser for help on using the repository browser.