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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5