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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (show annotations)
Thu Jan 3 13:16:27 1991 UTC (23 years, 3 months ago) by ram
Branch: MAIN
Changes since 1.5: +53 -37 lines
Fixed EXPT to correctly return complex results for negative numbers
raised to fractional powers.
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.6 1991/01/03 13:16:27 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 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
27 ;;; call, and saves number consing to boot.
28 ;;;
29 (defmacro def-math-rtn (name num-args)
30 (let ((function (intern (concatenate 'simple-string
31 "%"
32 (string-upcase name)))))
33 `(progn
34 (proclaim '(inline ,function))
35 (export ',function)
36 (def-c-routine (,name ,function) (double-float)
37 ,@(let ((results nil))
38 (dotimes (i num-args (nreverse results))
39 (push (list (intern (format nil "ARG-~D" i))
40 'double-float)
41 results)))))))
42
43 (eval-when (compile load eval)
44
45 (defun handle-reals (function var)
46 `((((foreach fixnum single-float bignum ratio))
47 (coerce (,function (coerce ,var 'double-float)) 'single-float))
48 ((double-float)
49 (,function ,var))))
50
51 ); eval-when (compile load eval)
52
53
54 ;;;; Stubs for the Unix math library.
55
56 ;;; Please refer to the Unix man pages for details about these routines.
57
58 ;;; Trigonometric.
59 (def-math-rtn "sin" 1)
60 (def-math-rtn "cos" 1)
61 (def-math-rtn "tan" 1)
62 (def-math-rtn "asin" 1)
63 (def-math-rtn "acos" 1)
64 (def-math-rtn "atan" 1)
65 (def-math-rtn "atan2" 2)
66 (def-math-rtn "sinh" 1)
67 (def-math-rtn "cosh" 1)
68 (def-math-rtn "tanh" 1)
69 (def-math-rtn "asinh" 1)
70 (def-math-rtn "acosh" 1)
71 (def-math-rtn "atanh" 1)
72
73 ;;; Exponential and Logarithmic.
74 (def-math-rtn "exp" 1)
75 (def-math-rtn "expm1" 1)
76 (def-math-rtn "log" 1)
77 (def-math-rtn "log10" 1)
78 (def-math-rtn "log1p" 1)
79 (def-math-rtn "pow" 2)
80 (def-math-rtn "cbrt" 1)
81 (def-math-rtn "sqrt" 1)
82 (def-math-rtn "hypot" 2)
83
84
85 ;;;; Power functions.
86
87 (defun exp (number)
88 "Return e raised to the power NUMBER."
89 (number-dispatch ((number number))
90 (handle-reals %exp number)
91 ((complex)
92 (* (exp (realpart number))
93 (cis (imagpart number))))))
94
95 ;;; INTEXP -- Handle the rational base, integer power case.
96
97 (defparameter *intexp-maximum-exponent* 10000)
98
99 ;;; This function precisely calculates base raised to an integral power. It
100 ;;; separates the cases by the sign of power, for efficiency reasons, as powers
101 ;;; can be calculated more efficiently if power is a positive integer. Values
102 ;;; of power are calculated as positive integers, and inverted if negative.
103 ;;;
104 (defun intexp (base power)
105 (when (> (abs power) *intexp-maximum-exponent*)
106 (cerror "Continue with calculation."
107 "The absolute value of ~S exceeds ~S."
108 power '*intexp-maximum-exponent* base power))
109 (cond ((minusp power)
110 (/ (intexp base (- power))))
111 ((eql base 2)
112 (ash 1 power))
113 (t
114 (do ((nextn (ash power -1) (ash power -1))
115 (total (if (oddp power) base 1)
116 (if (oddp power) (* base total) total)))
117 ((zerop nextn) total)
118 (setq base (* base base))
119 (setq power nextn)))))
120
121
122 ;;; EXPT -- Public
123 ;;;
124 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
125 ;;; floating point stuff. If both args are real, we try %POW right off,
126 ;;; assuming it will return 0 if the result may be complex. If so, we call
127 ;;; COMPLEX-POW which directly computes the complex result. We also separate
128 ;;; the complex-real and real-complex cases from the general complex case.
129 ;;;
130 (defun expt (base power)
131 "Returns BASE raised to the POWER."
132 (if (zerop power)
133 (1+ (* base power))
134 (labels ((real-expt (base power rtype)
135 (let* ((fbase (coerce base 'double-float))
136 (fpower (coerce power 'double-float))
137 (res (coerce (%pow fbase fpower) rtype)))
138 (if (and (zerop res) (minusp fbase))
139 (multiple-value-bind (re im)
140 (complex-pow fbase fpower)
141 (%make-complex (coerce re rtype) (coerce im rtype)))
142 res)))
143 (complex-pow (fbase fpower)
144 (let ((pow (%pow (- fbase) fpower))
145 (fpower*pi (* fpower pi)))
146 (values (* pow (%cos fpower*pi))
147 (* pow (%sin fpower*pi))))))
148 (declare (inline real-expt))
149 (number-dispatch ((base number) (power number))
150 (((foreach fixnum (or bignum ratio) (complex rational)) integer)
151 (intexp base power))
152 (((foreach single-float double-float) integer)
153 (real-expt base power '(dispatch-type base)))
154 (((foreach fixnum (or bignum ratio) single-float)
155 (foreach ratio single-float))
156 (real-expt base power 'single-float))
157 (((foreach fixnum (or bignum ratio) single-float double-float)
158 double-float)
159 (real-expt base power 'double-float))
160 ((double-float single-float)
161 (real-expt base power 'double-float))
162 (((foreach (complex rational) (complex float)) rational)
163 (* (expt (abs base) power)
164 (cis (* power (phase base)))))
165 (((foreach fixnum (or bignum ratio) single-float double-float)
166 complex)
167 (if (minusp base)
168 (/ (exp (* power (truly-the float (log (- base))))))
169 (exp (* power (truly-the float (log base))))))
170 (((foreach (complex float) (complex rational)) complex)
171 (exp (* power (log base))))))))
172
173 (defun log (number &optional (base nil base-p))
174 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
175 (if base-p
176 (/ (log number) (log base))
177 (number-dispatch ((number number))
178 (((foreach fixnum bignum ratio single-float))
179 (if (minusp number)
180 (complex (log (- number)) (coerce pi 'single-float))
181 (coerce (%log (coerce number 'double-float)) 'single-float)))
182 ((double-float)
183 (if (minusp number)
184 (complex (log (- number)) (coerce pi 'double-float))
185 (%log number)))
186 ((complex) (complex (log (abs number)) (phase number))))))
187
188 (defun sqrt (number)
189 "Return the square root of NUMBER."
190 (number-dispatch ((number number))
191 (((foreach fixnum bignum ratio single-float))
192 (if (minusp number)
193 (exp (/ (log number) 2))
194 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
195 ((double-float)
196 (if (minusp number)
197 (exp (/ (log number) 2))
198 (%sqrt number)))
199 ((complex) (exp (/ (log number) 2)))))
200
201 ;;; ISQRT: Integer square root - isqrt(n)**2 <= n
202 ;;; Upper and lower bounds on the result are estimated using integer-length.
203 ;;; On each iteration, one of the bounds is replaced by their mean.
204 ;;; The lower bound is returned when the bounds meet or differ by only 1.
205 ;;; Initial bounds guarantee that lg(sqrt(n)) = lg(n)/2 iterations suffice.
206
207 (defun isqrt (n)
208 "Returns the root of the nearest integer less than
209 n which is a perfect square."
210 (if (and (integerp n) (not (minusp n)))
211 (do* ((lg (integer-length n))
212 (lo (ash 1 (ash (1- lg) -1)))
213 (hi (+ lo (ash lo (if (oddp lg) -1 0))))) ;tighten by 3/4 if possible.
214 ((<= (1- hi) lo) lo)
215 (let ((mid (ash (+ lo hi) -1)))
216 (if (<= (* mid mid) n) (setq lo mid) (setq hi mid))))
217 (error "Isqrt: ~S argument must be a nonnegative integer" n)))
218
219
220
221 ;;;; Trigonometic and Related Functions
222
223 (defun abs (number)
224 "Returns the absolute value of the number."
225 (number-dispatch ((number number))
226 (((foreach single-float double-float fixnum rational))
227 (abs number))
228 ((complex)
229 (let ((rx (realpart number))
230 (ix (imagpart number)))
231 (etypecase rx
232 (rational
233 (sqrt (+ (* rx rx) (* ix ix))))
234 (single-float
235 (coerce (%hypot (coerce rx 'double-float)
236 (coerce ix 'double-float))
237 'single-float))
238 (double-float
239 (%hypot rx ix)))))))
240
241 (defun phase (number)
242 "Returns the angle part of the polar representation of a complex number.
243 For complex numbers, this is (atan (imagpart number) (realpart number)).
244 For non-complex positive numbers, this is 0. For non-complex negative
245 numbers this is PI."
246 (etypecase number
247 ((or rational single-float)
248 (if (minusp number)
249 (coerce pi 'single-float)
250 0.0f0))
251 (double-float
252 (if (minusp number)
253 (coerce pi 'double-float)
254 0.0d0))
255 (complex
256 (atan (imagpart number) (realpart number)))))
257
258
259 (defun sin (number)
260 "Return the sine of NUMBER."
261 (number-dispatch ((number number))
262 (handle-reals %sin number)
263 ((complex)
264 (let ((x (realpart number))
265 (y (imagpart number)))
266 (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))))
267
268 (defun cos (number)
269 "Return the cosine of NUMBER."
270 (number-dispatch ((number number))
271 (handle-reals %cos number)
272 ((complex)
273 (let ((x (realpart number))
274 (y (imagpart number)))
275 (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))))
276
277 (defun tan (number)
278 "Return the tangent of NUMBER."
279 (number-dispatch ((number number))
280 (handle-reals %tan number)
281 ((complex)
282 (let* ((num (sin number))
283 (denom (cos number)))
284 (if (zerop denom) (error "~S undefined tangent." number)
285 (/ num denom))))))
286
287 (defun cis (theta)
288 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
289 (if (complexp theta)
290 (error "Argument to CIS is complex: ~S" theta)
291 (complex (cos theta) (sin theta))))
292
293 (proclaim '(inline mult-by-i))
294 (defun mult-by-i (number)
295 (complex (imagpart number)
296 (- (realpart number))))
297
298 (defun complex-asin (number)
299 (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
300
301 (defun asin (number)
302 "Return the arc sine of NUMBER."
303 (number-dispatch ((number number))
304 ((rational)
305 (if (or (> number 1) (< number -1))
306 (complex-asin number)
307 (coerce (%asin (coerce number 'double-float)) 'single-float)))
308 (((foreach single-float double-float))
309 (if (or (> number (coerce 1 '(dispatch-type number)))
310 (< number (coerce -1 '(dispatch-type number))))
311 (complex-asin number)
312 (coerce (%asin (coerce number 'double-float))
313 '(dispatch-type number))))
314 ((complex)
315 (complex-asin number))))
316
317 (defun complex-acos (number)
318 (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))
319
320 (defun acos (number)
321 "Return the arc cosine of NUMBER."
322 (number-dispatch ((number number))
323 ((rational)
324 (if (or (> number 1) (< number -1))
325 (complex-acos number)
326 (coerce (%acos (coerce number 'double-float)) 'single-float)))
327 (((foreach single-float double-float))
328 (if (or (> number (coerce 1 '(dispatch-type number)))
329 (< number (coerce -1 '(dispatch-type number))))
330 (complex-acos number)
331 (coerce (%acos (coerce number 'double-float))
332 '(dispatch-type number))))
333 ((complex)
334 (complex-acos number))))
335
336
337 (defun atan (y &optional (x nil xp))
338 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
339 (if xp
340 (if (and (zerop x) (zerop y))
341 (if (plusp (float-sign (float x)))
342 y
343 (if (minusp (float-sign (float y)))
344 (- pi)
345 pi))
346 (number-dispatch ((y real) (x real))
347 (((foreach fixnum bignum ratio single-float)
348 (foreach fixnum bignum ratio single-float))
349 (coerce (%atan2 (coerce y 'double-float)
350 (coerce x 'double-float))
351 'single-float))
352 ((double-float (foreach fixnum bignum ratio single-float))
353 (%atan2 y (coerce x 'double-float)))
354 (((foreach fixnum bignum ratio single-float double-float)
355 double-float)
356 (%atan2 (coerce y 'double-float) x))))
357 (number-dispatch ((y number))
358 (handle-reals %atan y)
359 ((complex)
360 (let ((im (imagpart y))
361 (re (realpart y)))
362 (/ (- (log (complex (- 1 im) re))
363 (log (complex (+ 1 im) (- re))))
364 (complex 0 2)))))))
365
366
367 (defun sinh (number)
368 "Return the hyperbolic sine of NUMBER."
369 (/ (- (exp number) (exp (- number))) 2))
370
371 (defun cosh (number)
372 "Return the hyperbolic cosine of NUMBER."
373 (/ (+ (exp number) (exp (- number))) 2))
374
375 (defun tanh (number)
376 "Return the hyperbolic tangent of NUMBER."
377 (/ (- (exp number) (exp (- number)))
378 (+ (exp number) (exp (- number)))))
379
380
381 (defun asinh (number)
382 "Return the hyperbolic arc sine of NUMBER."
383 (log (+ number (sqrt (1+ (* number number))))))
384
385 (defun acosh (number)
386 "Return the hyperbolic arc cosine of NUMBER."
387 (log (+ number (* (1+ number) (sqrt (/ (1- number) (1+ number)))))))
388
389 (defun atanh (number)
390 "Return the hyperbolic arc tangent of NUMBER."
391 (log (* (1+ number) (sqrt (/ (- 1 (* number number)))))))
392

  ViewVC Help
Powered by ViewVC 1.1.5