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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.16 - (hide annotations)
Mon Jul 8 19:17:28 1996 UTC (17 years, 9 months ago) by ram
Branch: MAIN
Changes since 1.15: +6 -7 lines
Merge fix from Douglas Crosher to make acosh and atanh return the right thing.
1 wlott 1.1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2     ;;;
3     ;;; **********************************************************************
4 ram 1.7 ;;; 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     ;;;
7     (ext:file-comment
8 ram 1.16 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.16 1996/07/08 19:17:28 ram Exp $")
9 ram 1.7 ;;;
10 wlott 1.1 ;;; **********************************************************************
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 wlott 1.2 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
25 wlott 1.1
26 ram 1.5 ;;; 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 wlott 1.1 (defmacro def-math-rtn (name num-args)
30     (let ((function (intern (concatenate 'simple-string
31     "%"
32     (string-upcase name)))))
33 ram 1.4 `(progn
34 ram 1.5 (proclaim '(inline ,function))
35 ram 1.4 (export ',function)
36 wlott 1.10 (alien:def-alien-routine (,name ,function) double-float
37 ram 1.4 ,@(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 wlott 1.1
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 wlott 1.14 #-hpux (def-math-rtn "asinh" 1)
70     #-hpux (def-math-rtn "acosh" 1)
71     #-hpux (def-math-rtn "atanh" 1)
72 wlott 1.1
73     ;;; Exponential and Logarithmic.
74     (def-math-rtn "exp" 1)
75     (def-math-rtn "log" 1)
76     (def-math-rtn "log10" 1)
77     (def-math-rtn "pow" 2)
78     (def-math-rtn "sqrt" 1)
79     (def-math-rtn "hypot" 2)
80    
81    
82     ;;;; Power functions.
83    
84     (defun exp (number)
85     "Return e raised to the power NUMBER."
86     (number-dispatch ((number number))
87     (handle-reals %exp number)
88     ((complex)
89     (* (exp (realpart number))
90     (cis (imagpart number))))))
91    
92     ;;; INTEXP -- Handle the rational base, integer power case.
93    
94     (defparameter *intexp-maximum-exponent* 10000)
95    
96 ram 1.6 ;;; This function precisely calculates base raised to an integral power. It
97     ;;; separates the cases by the sign of power, for efficiency reasons, as powers
98     ;;; can be calculated more efficiently if power is a positive integer. Values
99     ;;; of power are calculated as positive integers, and inverted if negative.
100     ;;;
101 wlott 1.1 (defun intexp (base power)
102     (when (> (abs power) *intexp-maximum-exponent*)
103     (cerror "Continue with calculation."
104     "The absolute value of ~S exceeds ~S."
105     power '*intexp-maximum-exponent* base power))
106     (cond ((minusp power)
107     (/ (intexp base (- power))))
108     ((eql base 2)
109     (ash 1 power))
110     (t
111     (do ((nextn (ash power -1) (ash power -1))
112     (total (if (oddp power) base 1)
113     (if (oddp power) (* base total) total)))
114     ((zerop nextn) total)
115     (setq base (* base base))
116     (setq power nextn)))))
117    
118    
119 ram 1.6 ;;; EXPT -- Public
120     ;;;
121     ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
122     ;;; floating point stuff. If both args are real, we try %POW right off,
123     ;;; assuming it will return 0 if the result may be complex. If so, we call
124     ;;; COMPLEX-POW which directly computes the complex result. We also separate
125     ;;; the complex-real and real-complex cases from the general complex case.
126     ;;;
127 wlott 1.1 (defun expt (base power)
128 wlott 1.3 "Returns BASE raised to the POWER."
129 wlott 1.1 (if (zerop power)
130 ram 1.6 (1+ (* base power))
131     (labels ((real-expt (base power rtype)
132     (let* ((fbase (coerce base 'double-float))
133     (fpower (coerce power 'double-float))
134     (res (coerce (%pow fbase fpower) rtype)))
135     (if (and (zerop res) (minusp fbase))
136     (multiple-value-bind (re im)
137     (complex-pow fbase fpower)
138     (%make-complex (coerce re rtype) (coerce im rtype)))
139     res)))
140     (complex-pow (fbase fpower)
141     (let ((pow (%pow (- fbase) fpower))
142     (fpower*pi (* fpower pi)))
143     (values (* pow (%cos fpower*pi))
144     (* pow (%sin fpower*pi))))))
145     (declare (inline real-expt))
146     (number-dispatch ((base number) (power number))
147     (((foreach fixnum (or bignum ratio) (complex rational)) integer)
148     (intexp base power))
149 ram 1.8 (((foreach single-float double-float) rational)
150 ram 1.6 (real-expt base power '(dispatch-type base)))
151 ram 1.9 (((foreach fixnum (or bignum ratio) single-float)
152 ram 1.6 (foreach ratio single-float))
153     (real-expt base power 'single-float))
154     (((foreach fixnum (or bignum ratio) single-float double-float)
155     double-float)
156     (real-expt base power 'double-float))
157     ((double-float single-float)
158     (real-expt base power 'double-float))
159     (((foreach (complex rational) (complex float)) rational)
160     (* (expt (abs base) power)
161     (cis (* power (phase base)))))
162     (((foreach fixnum (or bignum ratio) single-float double-float)
163     complex)
164     (if (minusp base)
165     (/ (exp (* power (truly-the float (log (- base))))))
166     (exp (* power (truly-the float (log base))))))
167     (((foreach (complex float) (complex rational)) complex)
168     (exp (* power (log base))))))))
169 wlott 1.1
170     (defun log (number &optional (base nil base-p))
171 wlott 1.3 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
172 wlott 1.1 (if base-p
173     (/ (log number) (log base))
174     (number-dispatch ((number number))
175 wlott 1.3 (((foreach fixnum bignum ratio single-float))
176     (if (minusp number)
177     (complex (log (- number)) (coerce pi 'single-float))
178     (coerce (%log (coerce number 'double-float)) 'single-float)))
179     ((double-float)
180     (if (minusp number)
181     (complex (log (- number)) (coerce pi 'double-float))
182     (%log number)))
183     ((complex) (complex (log (abs number)) (phase number))))))
184 wlott 1.1
185     (defun sqrt (number)
186     "Return the square root of NUMBER."
187     (number-dispatch ((number number))
188 wlott 1.3 (((foreach fixnum bignum ratio single-float))
189     (if (minusp number)
190     (exp (/ (log number) 2))
191     (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
192     ((double-float)
193     (if (minusp number)
194     (exp (/ (log number) 2))
195     (%sqrt number)))
196     ((complex) (exp (/ (log number) 2)))))
197 wlott 1.1
198    
199     ;;;; Trigonometic and Related Functions
200    
201 wlott 1.2 (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 wlott 1.1
219     (defun phase (number)
220     "Returns the angle part of the polar representation of a complex number.
221 wlott 1.3 For complex numbers, this is (atan (imagpart number) (realpart number)).
222 wlott 1.1 For non-complex positive numbers, this is 0. For non-complex negative
223     numbers this is PI."
224 wlott 1.3 (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 wlott 1.1
236    
237     (defun sin (number)
238 wlott 1.3 "Return the sine of NUMBER."
239 wlott 1.1 (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 wlott 1.3 "Return the cosine of NUMBER."
248 wlott 1.1 (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 wlott 1.3 "Return the tangent of NUMBER."
257 wlott 1.1 (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 wlott 1.2 (defun cis (theta)
266 wlott 1.3 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
267 wlott 1.2 (if (complexp theta)
268     (error "Argument to CIS is complex: ~S" theta)
269     (complex (cos theta) (sin theta))))
270 wlott 1.1
271 wlott 1.3 (proclaim '(inline mult-by-i))
272     (defun mult-by-i (number)
273     (complex (imagpart number)
274     (- (realpart number))))
275 wlott 1.1
276 wlott 1.3 (defun complex-asin (number)
277     (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
278    
279 wlott 1.1 (defun asin (number)
280 wlott 1.3 "Return the arc sine of NUMBER."
281 wlott 1.1 (number-dispatch ((number number))
282 wlott 1.3 ((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 wlott 1.1 ((complex)
293 wlott 1.3 (complex-asin number))))
294 wlott 1.1
295 wlott 1.3 (defun complex-acos (number)
296     (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))
297 wlott 1.1
298     (defun acos (number)
299 wlott 1.3 "Return the arc cosine of NUMBER."
300 wlott 1.1 (number-dispatch ((number number))
301 wlott 1.3 ((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 wlott 1.1 ((complex)
312 wlott 1.3 (complex-acos number))))
313 wlott 1.1
314    
315     (defun atan (y &optional (x nil xp))
316 wlott 1.3 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
317 wlott 1.1 (if xp
318 wlott 1.12 (flet ((atan2 (y x)
319 wlott 1.13 (declare (type double-float y x)
320     (values double-float))
321     (if (zerop x)
322     (if (zerop y)
323     (if (plusp (float-sign x))
324     y
325     (float-sign y pi))
326     (float-sign y (/ pi 2)))
327     (%atan2 y x))))
328     (number-dispatch ((y number) (x number))
329     ((double-float
330     (foreach double-float single-float fixnum bignum ratio))
331     (atan2 y (coerce x 'double-float)))
332     (((foreach single-float fixnum bignum ratio)
333     double-float)
334     (atan2 (coerce y 'double-float) x))
335     (((foreach single-float fixnum bignum ratio)
336     (foreach single-float fixnum bignum ratio))
337     (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
338     'single-float))))
339 wlott 1.1 (number-dispatch ((y number))
340     (handle-reals %atan y)
341     ((complex)
342 wlott 1.3 (let ((im (imagpart y))
343     (re (realpart y)))
344     (/ (- (log (complex (- 1 im) re))
345     (log (complex (+ 1 im) (- re))))
346     (complex 0 2)))))))
347    
348    
349     (defun sinh (number)
350     "Return the hyperbolic sine of NUMBER."
351     (/ (- (exp number) (exp (- number))) 2))
352    
353     (defun cosh (number)
354     "Return the hyperbolic cosine of NUMBER."
355     (/ (+ (exp number) (exp (- number))) 2))
356    
357     (defun tanh (number)
358     "Return the hyperbolic tangent of NUMBER."
359     (/ (- (exp number) (exp (- number)))
360     (+ (exp number) (exp (- number)))))
361    
362    
363     (defun asinh (number)
364     "Return the hyperbolic arc sine of NUMBER."
365     (log (+ number (sqrt (1+ (* number number))))))
366    
367     (defun acosh (number)
368     "Return the hyperbolic arc cosine of NUMBER."
369 ram 1.16 (* 2 (log (+ (sqrt (/ (1+ number) 2)) (sqrt (/ (1- number) 2))))))
370 wlott 1.3
371     (defun atanh (number)
372     "Return the hyperbolic arc tangent of NUMBER."
373 ram 1.16 (/ (- (log (1+ number)) (log (- 1 number))) 2))
374 wlott 1.3
375 wlott 1.14
376     ;;; HP-UX does not supply C versions of asinh, acosh, and atanh, so just
377     ;;; use the same equations as above, but force double-float arith.
378    
379     #+hpux
380     (defun %asinh (number)
381     (declare (type double-float number)
382     (values double-float)
383     (optimize (speed 3) (safety 0) (inhibit-warnings 3)))
384     (log (+ number (the double-float (sqrt (1+ (* number number)))))))
385    
386     #+hpux
387     (defun %acosh (number)
388     (declare (type double-float number)
389     (values double-float)
390     (optimize (speed 3) (safety 0) (inhibit-warnings 3)))
391 ram 1.16 (* 2 (log (+ (the double-float (sqrt (/ (1+ number) 2)))
392     (the double-float (sqrt (/ (1- number) 2)))))))
393 wlott 1.14
394     #+hpux
395     (defun %atanh (number)
396     (declare (type double-float number)
397     (values double-float)
398     (optimize (speed 3) (safety 0) (inhibit-warnings 3)))
399 ram 1.16 (/ (- (log (1+ number)) (log (- 1 number))) 2))

  ViewVC Help
Powered by ViewVC 1.1.5