/[cmucl]/src/code/irrat.lisp
ViewVC logotype

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (show annotations)
Thu Oct 18 02:56:15 1990 UTC (23 years, 6 months ago) by ram
Branch: MAIN
Changes since 1.3: +20 -20 lines
Made all the C stubs be exported and maybe-inline.  Put the generic
functions under (space 0) so that the stubs will be inline expanded.
Changed some uses of INTEGER-DECODE-FLOAT to FLOAT-SIGN, since only
the sign was being used.
Note that EXPT is broken, since (expt -2.0 0.5) should be a complex,
not 0.
1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2 ;;;
3 ;;; **********************************************************************
4 ;;; This code was written as part of the Spice Lisp project at
5 ;;; Carnegie-Mellon University, and has been placed in the public domain.
6 ;;; If you want to use this code or any part of Spice Lisp, please contact
7 ;;; Scott Fahlman (FAHLMAN@CMUC).
8 ;;; **********************************************************************
9 ;;;
10 ;;; $Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.4 1990/10/18 02:56:15 ram Exp $
11 ;;;
12 ;;; This file contains all the irrational functions. Actually, most of the
13 ;;; work is done by calling out to C...
14 ;;;
15 ;;; Author: William Lott.
16 ;;;
17
18 (in-package "KERNEL")
19
20
21 ;;;; Random constants, utility functions, and macros.
22
23 (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
24 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
25
26 (defmacro def-math-rtn (name num-args)
27 (let ((function (intern (concatenate 'simple-string
28 "%"
29 (string-upcase name)))))
30 `(progn
31 (proclaim '(maybe-inline ,function))
32 (export ',function)
33 (def-c-routine (,name ,function) (double-float)
34 ,@(let ((results nil))
35 (dotimes (i num-args (nreverse results))
36 (push (list (intern (format nil "ARG-~D" i))
37 'double-float)
38 results)))))))
39
40 (eval-when (compile load eval)
41
42 (defun handle-reals (function var)
43 `((((foreach fixnum single-float bignum ratio))
44 (coerce (,function (coerce ,var 'double-float)) 'single-float))
45 ((double-float)
46 (,function ,var))))
47
48 ); eval-when (compile load eval)
49
50
51 ;;;; Stubs for the Unix math library.
52
53 ;;; Please refer to the Unix man pages for details about these routines.
54
55 ;;; Trigonometric.
56 (def-math-rtn "sin" 1)
57 (def-math-rtn "cos" 1)
58 (def-math-rtn "tan" 1)
59 (def-math-rtn "asin" 1)
60 (def-math-rtn "acos" 1)
61 (def-math-rtn "atan" 1)
62 (def-math-rtn "atan2" 2)
63 (def-math-rtn "sinh" 1)
64 (def-math-rtn "cosh" 1)
65 (def-math-rtn "tanh" 1)
66 (def-math-rtn "asinh" 1)
67 (def-math-rtn "acosh" 1)
68 (def-math-rtn "atanh" 1)
69
70 ;;; Exponential and Logarithmic.
71 (def-math-rtn "exp" 1)
72 (def-math-rtn "expm1" 1)
73 (def-math-rtn "log" 1)
74 (def-math-rtn "log10" 1)
75 (def-math-rtn "log1p" 1)
76 (def-math-rtn "pow" 2)
77 (def-math-rtn "cbrt" 1)
78 (def-math-rtn "sqrt" 1)
79 (def-math-rtn "hypot" 2)
80
81
82 ;;;; Power functions.
83
84 ;;; Let the C stubs be opportunistically inline expanded.
85 ;;;
86 (proclaim '(optimize (space 0)))
87
88 (defun exp (number)
89 "Return e raised to the power NUMBER."
90 (number-dispatch ((number number))
91 (handle-reals %exp number)
92 ((complex)
93 (* (exp (realpart number))
94 (cis (imagpart number))))))
95
96 ;;; INTEXP -- Handle the rational base, integer power case.
97
98 (defparameter *intexp-maximum-exponent* 10000)
99
100 (defun intexp (base power)
101 (when (> (abs power) *intexp-maximum-exponent*)
102 (cerror "Continue with calculation."
103 "The absolute value of ~S exceeds ~S."
104 power '*intexp-maximum-exponent* base power))
105 (cond ((minusp power)
106 (/ (intexp base (- power))))
107 ((eql base 2)
108 (ash 1 power))
109 (t
110 (do ((nextn (ash power -1) (ash power -1))
111 (total (if (oddp power) base 1)
112 (if (oddp power) (* base total) total)))
113 ((zerop nextn) total)
114 (setq base (* base base))
115 (setq power nextn)))))
116
117 ;;; This function calculates x raised to the nth power. It separates
118 ;;; the cases by the type of n, for efficiency reasons, as powers can
119 ;;; be calculated more efficiently if n is a positive integer, Therefore,
120 ;;; All integer values of n are calculated as positive integers, and
121 ;;; inverted if negative.
122
123 (defun expt (base power)
124 "Returns BASE raised to the POWER."
125 (if (zerop power)
126 ;; This is wrong if power isn't an integer.
127 (typecase (realpart base)
128 (single-float (coerce 1 'single-float))
129 (double-float (coerce 1 'double-float))
130 (t 1))
131 (number-dispatch ((base number) (power number))
132 (((foreach fixnum bignum ratio (complex rational)) integer)
133 (intexp base power))
134 (((foreach single-float double-float) integer)
135 (coerce (%pow (coerce base 'double-float)
136 (coerce power 'double-float))
137 '(dispatch-type base)))
138 (((foreach fixnum bignum ratio single-float)
139 (foreach ratio single-float))
140 (coerce (%pow (coerce base 'double-float)
141 (coerce power 'double-float))
142 'single-float))
143 (((foreach fixnum bignum ratio single-float double-float) double-float)
144 (%pow (coerce base 'double-float) (coerce power 'double-float)))
145 (((complex rational) ratio)
146 (* (expt (abs base) power)
147 (cis (* power (phase base)))))
148 (((complex float) (foreach integer ratio))
149 (* (expt (abs base) power)
150 (cis (* power (phase base)))))
151 (((foreach fixnum bignum ratio single-float double-float) complex)
152 (if (minusp base)
153 (/ (exp (* power (log (- base)))))
154 (exp (* power (log base)))))
155 (((foreach (complex float) (complex rational)) complex)
156 (exp (* power (log base)))))))
157
158 (defun log (number &optional (base nil base-p))
159 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
160 (if base-p
161 (/ (log number) (log base))
162 (number-dispatch ((number number))
163 (((foreach fixnum bignum ratio single-float))
164 (if (minusp number)
165 (complex (log (- number)) (coerce pi 'single-float))
166 (coerce (%log (coerce number 'double-float)) 'single-float)))
167 ((double-float)
168 (if (minusp number)
169 (complex (log (- number)) (coerce pi 'double-float))
170 (%log number)))
171 ((complex) (complex (log (abs number)) (phase number))))))
172
173 (defun sqrt (number)
174 "Return the square root of NUMBER."
175 (number-dispatch ((number number))
176 (((foreach fixnum bignum ratio single-float))
177 (if (minusp number)
178 (exp (/ (log number) 2))
179 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
180 ((double-float)
181 (if (minusp number)
182 (exp (/ (log number) 2))
183 (%sqrt number)))
184 ((complex) (exp (/ (log number) 2)))))
185
186 ;;; ISQRT: Integer square root - isqrt(n)**2 <= n
187 ;;; Upper and lower bounds on the result are estimated using integer-length.
188 ;;; On each iteration, one of the bounds is replaced by their mean.
189 ;;; The lower bound is returned when the bounds meet or differ by only 1.
190 ;;; Initial bounds guarantee that lg(sqrt(n)) = lg(n)/2 iterations suffice.
191
192 (defun isqrt (n)
193 "Returns the root of the nearest integer less than
194 n which is a perfect square."
195 (if (and (integerp n) (not (minusp n)))
196 (do* ((lg (integer-length n))
197 (lo (ash 1 (ash (1- lg) -1)))
198 (hi (+ lo (ash lo (if (oddp lg) -1 0))))) ;tighten by 3/4 if possible.
199 ((<= (1- hi) lo) lo)
200 (let ((mid (ash (+ lo hi) -1)))
201 (if (<= (* mid mid) n) (setq lo mid) (setq hi mid))))
202 (error "Isqrt: ~S argument must be a nonnegative integer" n)))
203
204
205
206 ;;;; Trigonometic and Related Functions
207
208 (defun abs (number)
209 "Returns the absolute value of the number."
210 (number-dispatch ((number number))
211 (((foreach single-float double-float fixnum rational))
212 (abs number))
213 ((complex)
214 (let ((rx (realpart number))
215 (ix (imagpart number)))
216 (etypecase rx
217 (rational
218 (sqrt (+ (* rx rx) (* ix ix))))
219 (single-float
220 (coerce (%hypot (coerce rx 'double-float)
221 (coerce ix 'double-float))
222 'single-float))
223 (double-float
224 (%hypot rx ix)))))))
225
226 (defun phase (number)
227 "Returns the angle part of the polar representation of a complex number.
228 For complex numbers, this is (atan (imagpart number) (realpart number)).
229 For non-complex positive numbers, this is 0. For non-complex negative
230 numbers this is PI."
231 (etypecase number
232 ((or rational single-float)
233 (if (minusp number)
234 (coerce pi 'single-float)
235 0.0f0))
236 (double-float
237 (if (minusp number)
238 (coerce pi 'double-float)
239 0.0d0))
240 (complex
241 (atan (imagpart number) (realpart number)))))
242
243
244 (defun sin (number)
245 "Return the sine of NUMBER."
246 (number-dispatch ((number number))
247 (handle-reals %sin number)
248 ((complex)
249 (let ((x (realpart number))
250 (y (imagpart number)))
251 (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))))
252
253 (defun cos (number)
254 "Return the cosine of NUMBER."
255 (number-dispatch ((number number))
256 (handle-reals %cos number)
257 ((complex)
258 (let ((x (realpart number))
259 (y (imagpart number)))
260 (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))))
261
262 (defun tan (number)
263 "Return the tangent of NUMBER."
264 (number-dispatch ((number number))
265 (handle-reals %tan number)
266 ((complex)
267 (let* ((num (sin number))
268 (denom (cos number)))
269 (if (zerop denom) (error "~S undefined tangent." number)
270 (/ num denom))))))
271
272 (defun cis (theta)
273 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
274 (if (complexp theta)
275 (error "Argument to CIS is complex: ~S" theta)
276 (complex (cos theta) (sin theta))))
277
278 (proclaim '(inline mult-by-i))
279 (defun mult-by-i (number)
280 (complex (imagpart number)
281 (- (realpart number))))
282
283 (defun complex-asin (number)
284 (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
285
286 (defun asin (number)
287 "Return the arc sine of NUMBER."
288 (number-dispatch ((number number))
289 ((rational)
290 (if (or (> number 1) (< number -1))
291 (complex-asin number)
292 (coerce (%asin (coerce number 'double-float)) 'single-float)))
293 (((foreach single-float double-float))
294 (if (or (> number (coerce 1 '(dispatch-type number)))
295 (< number (coerce -1 '(dispatch-type number))))
296 (complex-asin number)
297 (coerce (%asin (coerce number 'double-float))
298 '(dispatch-type number))))
299 ((complex)
300 (complex-asin number))))
301
302 (defun complex-acos (number)
303 (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))
304
305 (defun acos (number)
306 "Return the arc cosine of NUMBER."
307 (number-dispatch ((number number))
308 ((rational)
309 (if (or (> number 1) (< number -1))
310 (complex-acos number)
311 (coerce (%acos (coerce number 'double-float)) 'single-float)))
312 (((foreach single-float double-float))
313 (if (or (> number (coerce 1 '(dispatch-type number)))
314 (< number (coerce -1 '(dispatch-type number))))
315 (complex-acos number)
316 (coerce (%acos (coerce number 'double-float))
317 '(dispatch-type number))))
318 ((complex)
319 (complex-acos number))))
320
321
322 (defun atan (y &optional (x nil xp))
323 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
324 (if xp
325 (if (and (zerop x) (zerop y))
326 (if (plusp (float-sign (float x)))
327 y
328 (if (minusp (float-sign (float y)))
329 (- pi)
330 pi))
331 (number-dispatch ((y real) (x real))
332 (((foreach fixnum bignum ratio single-float)
333 (foreach fixnum bignum ratio single-float))
334 (coerce (%atan2 (coerce y 'double-float)
335 (coerce x 'double-float))
336 'single-float))
337 ((double-float (foreach fixnum bignum ratio single-float))
338 (%atan2 y (coerce x 'double-float)))
339 (((foreach fixnum bignum ratio single-float double-float)
340 double-float)
341 (%atan2 (coerce y 'double-float) x))))
342 (number-dispatch ((y number))
343 (handle-reals %atan y)
344 ((complex)
345 (let ((im (imagpart y))
346 (re (realpart y)))
347 (/ (- (log (complex (- 1 im) re))
348 (log (complex (+ 1 im) (- re))))
349 (complex 0 2)))))))
350
351
352 (defun sinh (number)
353 "Return the hyperbolic sine of NUMBER."
354 (/ (- (exp number) (exp (- number))) 2))
355
356 (defun cosh (number)
357 "Return the hyperbolic cosine of NUMBER."
358 (/ (+ (exp number) (exp (- number))) 2))
359
360 (defun tanh (number)
361 "Return the hyperbolic tangent of NUMBER."
362 (/ (- (exp number) (exp (- number)))
363 (+ (exp number) (exp (- number)))))
364
365
366 (defun asinh (number)
367 "Return the hyperbolic arc sine of NUMBER."
368 (log (+ number (sqrt (1+ (* number number))))))
369
370 (defun acosh (number)
371 "Return the hyperbolic arc cosine of NUMBER."
372 (log (+ number (* (1+ number) (sqrt (/ (1- number) (1+ number)))))))
373
374 (defun atanh (number)
375 "Return the hyperbolic arc tangent of NUMBER."
376 (log (* (1+ number) (sqrt (/ (- 1 (* number number)))))))
377
378 (proclaim '(optimize (space 1)))

  ViewVC Help
Powered by ViewVC 1.1.5