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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5