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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.11 - (show annotations)
Mon Mar 1 15:24:54 1993 UTC (21 years, 1 month ago) by ram
Branch: MAIN
Changes since 1.10: +1 -20 lines
Move ISQRT to numbers.lisp
1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2 ;;;
3 ;;; **********************************************************************
4 ;;; This code was written as part of the CMU Common 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 CMU Common Lisp, please contact
7 ;;; Scott Fahlman or slisp-group@cs.cmu.edu.
8 ;;;
9 (ext:file-comment
10 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.11 1993/03/01 15:24:54 ram Exp $")
11 ;;;
12 ;;; **********************************************************************
13 ;;;
14 ;;; This file contains all the irrational functions. Actually, most of the
15 ;;; work is done by calling out to C...
16 ;;;
17 ;;; Author: William Lott.
18 ;;;
19
20 (in-package "KERNEL")
21
22
23 ;;;; Random constants, utility functions, and macros.
24
25 (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
26 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
27
28 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
29 ;;; call, and saves number consing to boot.
30 ;;;
31 (defmacro def-math-rtn (name num-args)
32 (let ((function (intern (concatenate 'simple-string
33 "%"
34 (string-upcase name)))))
35 `(progn
36 (proclaim '(inline ,function))
37 (export ',function)
38 (alien:def-alien-routine (,name ,function) double-float
39 ,@(let ((results nil))
40 (dotimes (i num-args (nreverse results))
41 (push (list (intern (format nil "ARG-~D" i))
42 'double-float)
43 results)))))))
44
45 (eval-when (compile load eval)
46
47 (defun handle-reals (function var)
48 `((((foreach fixnum single-float bignum ratio))
49 (coerce (,function (coerce ,var 'double-float)) 'single-float))
50 ((double-float)
51 (,function ,var))))
52
53 ); eval-when (compile load eval)
54
55
56 ;;;; Stubs for the Unix math library.
57
58 ;;; Please refer to the Unix man pages for details about these routines.
59
60 ;;; Trigonometric.
61 (def-math-rtn "sin" 1)
62 (def-math-rtn "cos" 1)
63 (def-math-rtn "tan" 1)
64 (def-math-rtn "asin" 1)
65 (def-math-rtn "acos" 1)
66 (def-math-rtn "atan" 1)
67 (def-math-rtn "atan2" 2)
68 (def-math-rtn "sinh" 1)
69 (def-math-rtn "cosh" 1)
70 (def-math-rtn "tanh" 1)
71 (def-math-rtn "asinh" 1)
72 (def-math-rtn "acosh" 1)
73 (def-math-rtn "atanh" 1)
74
75 ;;; Exponential and Logarithmic.
76 (def-math-rtn "exp" 1)
77 (def-math-rtn "expm1" 1)
78 (def-math-rtn "log" 1)
79 (def-math-rtn "log10" 1)
80 (def-math-rtn "log1p" 1)
81 (def-math-rtn "pow" 2)
82 (def-math-rtn "cbrt" 1)
83 (def-math-rtn "sqrt" 1)
84 (def-math-rtn "hypot" 2)
85
86
87 ;;;; Power functions.
88
89 (defun exp (number)
90 "Return e raised to the power NUMBER."
91 (number-dispatch ((number number))
92 (handle-reals %exp number)
93 ((complex)
94 (* (exp (realpart number))
95 (cis (imagpart number))))))
96
97 ;;; INTEXP -- Handle the rational base, integer power case.
98
99 (defparameter *intexp-maximum-exponent* 10000)
100
101 ;;; This function precisely calculates base raised to an integral power. It
102 ;;; separates the cases by the sign of power, for efficiency reasons, as powers
103 ;;; can be calculated more efficiently if power is a positive integer. Values
104 ;;; of power are calculated as positive integers, and inverted if negative.
105 ;;;
106 (defun intexp (base power)
107 (when (> (abs power) *intexp-maximum-exponent*)
108 (cerror "Continue with calculation."
109 "The absolute value of ~S exceeds ~S."
110 power '*intexp-maximum-exponent* base power))
111 (cond ((minusp power)
112 (/ (intexp base (- power))))
113 ((eql base 2)
114 (ash 1 power))
115 (t
116 (do ((nextn (ash power -1) (ash power -1))
117 (total (if (oddp power) base 1)
118 (if (oddp power) (* base total) total)))
119 ((zerop nextn) total)
120 (setq base (* base base))
121 (setq power nextn)))))
122
123
124 ;;; EXPT -- Public
125 ;;;
126 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
127 ;;; floating point stuff. If both args are real, we try %POW right off,
128 ;;; assuming it will return 0 if the result may be complex. If so, we call
129 ;;; COMPLEX-POW which directly computes the complex result. We also separate
130 ;;; the complex-real and real-complex cases from the general complex case.
131 ;;;
132 (defun expt (base power)
133 "Returns BASE raised to the POWER."
134 (if (zerop power)
135 (1+ (* base power))
136 (labels ((real-expt (base power rtype)
137 (let* ((fbase (coerce base 'double-float))
138 (fpower (coerce power 'double-float))
139 (res (coerce (%pow fbase fpower) rtype)))
140 (if (and (zerop res) (minusp fbase))
141 (multiple-value-bind (re im)
142 (complex-pow fbase fpower)
143 (%make-complex (coerce re rtype) (coerce im rtype)))
144 res)))
145 (complex-pow (fbase fpower)
146 (let ((pow (%pow (- fbase) fpower))
147 (fpower*pi (* fpower pi)))
148 (values (* pow (%cos fpower*pi))
149 (* pow (%sin fpower*pi))))))
150 (declare (inline real-expt))
151 (number-dispatch ((base number) (power number))
152 (((foreach fixnum (or bignum ratio) (complex rational)) integer)
153 (intexp base power))
154 (((foreach single-float double-float) rational)
155 (real-expt base power '(dispatch-type base)))
156 (((foreach fixnum (or bignum ratio) single-float)
157 (foreach ratio single-float))
158 (real-expt base power 'single-float))
159 (((foreach fixnum (or bignum ratio) single-float double-float)
160 double-float)
161 (real-expt base power 'double-float))
162 ((double-float single-float)
163 (real-expt base power 'double-float))
164 (((foreach (complex rational) (complex float)) rational)
165 (* (expt (abs base) power)
166 (cis (* power (phase base)))))
167 (((foreach fixnum (or bignum ratio) single-float double-float)
168 complex)
169 (if (minusp base)
170 (/ (exp (* power (truly-the float (log (- base))))))
171 (exp (* power (truly-the float (log base))))))
172 (((foreach (complex float) (complex rational)) complex)
173 (exp (* power (log base))))))))
174
175 (defun log (number &optional (base nil base-p))
176 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
177 (if base-p
178 (/ (log number) (log base))
179 (number-dispatch ((number number))
180 (((foreach fixnum bignum ratio single-float))
181 (if (minusp number)
182 (complex (log (- number)) (coerce pi 'single-float))
183 (coerce (%log (coerce number 'double-float)) 'single-float)))
184 ((double-float)
185 (if (minusp number)
186 (complex (log (- number)) (coerce pi 'double-float))
187 (%log number)))
188 ((complex) (complex (log (abs number)) (phase number))))))
189
190 (defun sqrt (number)
191 "Return the square root of NUMBER."
192 (number-dispatch ((number number))
193 (((foreach fixnum bignum ratio single-float))
194 (if (minusp number)
195 (exp (/ (log number) 2))
196 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
197 ((double-float)
198 (if (minusp number)
199 (exp (/ (log number) 2))
200 (%sqrt number)))
201 ((complex) (exp (/ (log number) 2)))))
202
203
204 ;;;; Trigonometic and Related Functions
205
206 (defun abs (number)
207 "Returns the absolute value of the number."
208 (number-dispatch ((number number))
209 (((foreach single-float double-float fixnum rational))
210 (abs number))
211 ((complex)
212 (let ((rx (realpart number))
213 (ix (imagpart number)))
214 (etypecase rx
215 (rational
216 (sqrt (+ (* rx rx) (* ix ix))))
217 (single-float
218 (coerce (%hypot (coerce rx 'double-float)
219 (coerce ix 'double-float))
220 'single-float))
221 (double-float
222 (%hypot rx ix)))))))
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 ((or rational single-float)
231 (if (minusp number)
232 (coerce pi 'single-float)
233 0.0f0))
234 (double-float
235 (if (minusp number)
236 (coerce pi 'double-float)
237 0.0d0))
238 (complex
239 (atan (imagpart number) (realpart number)))))
240
241
242 (defun sin (number)
243 "Return the sine of NUMBER."
244 (number-dispatch ((number number))
245 (handle-reals %sin number)
246 ((complex)
247 (let ((x (realpart number))
248 (y (imagpart number)))
249 (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))))
250
251 (defun cos (number)
252 "Return the cosine of NUMBER."
253 (number-dispatch ((number number))
254 (handle-reals %cos number)
255 ((complex)
256 (let ((x (realpart number))
257 (y (imagpart number)))
258 (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))))
259
260 (defun tan (number)
261 "Return the tangent of NUMBER."
262 (number-dispatch ((number number))
263 (handle-reals %tan number)
264 ((complex)
265 (let* ((num (sin number))
266 (denom (cos number)))
267 (if (zerop denom) (error "~S undefined tangent." number)
268 (/ num denom))))))
269
270 (defun cis (theta)
271 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
272 (if (complexp theta)
273 (error "Argument to CIS is complex: ~S" theta)
274 (complex (cos theta) (sin theta))))
275
276 (proclaim '(inline mult-by-i))
277 (defun mult-by-i (number)
278 (complex (imagpart number)
279 (- (realpart number))))
280
281 (defun complex-asin (number)
282 (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
283
284 (defun asin (number)
285 "Return the arc sine of NUMBER."
286 (number-dispatch ((number number))
287 ((rational)
288 (if (or (> number 1) (< number -1))
289 (complex-asin number)
290 (coerce (%asin (coerce number 'double-float)) 'single-float)))
291 (((foreach single-float double-float))
292 (if (or (> number (coerce 1 '(dispatch-type number)))
293 (< number (coerce -1 '(dispatch-type number))))
294 (complex-asin number)
295 (coerce (%asin (coerce number 'double-float))
296 '(dispatch-type number))))
297 ((complex)
298 (complex-asin number))))
299
300 (defun complex-acos (number)
301 (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))
302
303 (defun acos (number)
304 "Return the arc cosine of NUMBER."
305 (number-dispatch ((number number))
306 ((rational)
307 (if (or (> number 1) (< number -1))
308 (complex-acos number)
309 (coerce (%acos (coerce number 'double-float)) 'single-float)))
310 (((foreach single-float double-float))
311 (if (or (> number (coerce 1 '(dispatch-type number)))
312 (< number (coerce -1 '(dispatch-type number))))
313 (complex-acos number)
314 (coerce (%acos (coerce number 'double-float))
315 '(dispatch-type number))))
316 ((complex)
317 (complex-acos number))))
318
319
320 (defun atan (y &optional (x nil xp))
321 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
322 (if xp
323 (if (and (zerop x) (zerop y))
324 (if (plusp (float-sign (float x)))
325 y
326 (if (minusp (float-sign (float y)))
327 (- pi)
328 pi))
329 (number-dispatch ((y real) (x real))
330 (((foreach fixnum bignum ratio single-float)
331 (foreach fixnum bignum ratio single-float))
332 (coerce (%atan2 (coerce y 'double-float)
333 (coerce x 'double-float))
334 'single-float))
335 ((double-float (foreach fixnum bignum ratio single-float))
336 (%atan2 y (coerce x 'double-float)))
337 (((foreach fixnum bignum ratio single-float double-float)
338 double-float)
339 (%atan2 (coerce y 'double-float) x))))
340 (number-dispatch ((y number))
341 (handle-reals %atan y)
342 ((complex)
343 (let ((im (imagpart y))
344 (re (realpart y)))
345 (/ (- (log (complex (- 1 im) re))
346 (log (complex (+ 1 im) (- re))))
347 (complex 0 2)))))))
348
349
350 (defun sinh (number)
351 "Return the hyperbolic sine of NUMBER."
352 (/ (- (exp number) (exp (- number))) 2))
353
354 (defun cosh (number)
355 "Return the hyperbolic cosine of NUMBER."
356 (/ (+ (exp number) (exp (- number))) 2))
357
358 (defun tanh (number)
359 "Return the hyperbolic tangent of NUMBER."
360 (/ (- (exp number) (exp (- number)))
361 (+ (exp number) (exp (- number)))))
362
363
364 (defun asinh (number)
365 "Return the hyperbolic arc sine of NUMBER."
366 (log (+ number (sqrt (1+ (* number number))))))
367
368 (defun acosh (number)
369 "Return the hyperbolic arc cosine of NUMBER."
370 (log (+ number (* (1+ number) (sqrt (/ (1- number) (1+ number)))))))
371
372 (defun atanh (number)
373 "Return the hyperbolic arc tangent of NUMBER."
374 (log (* (1+ number) (sqrt (/ (- 1 (* number number)))))))
375

  ViewVC Help
Powered by ViewVC 1.1.5