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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.12 - (hide annotations)
Sat May 8 04:49:36 1993 UTC (20 years, 11 months ago) by wlott
Branch: MAIN
Changes since 1.11: +36 -18 lines
Fixed atan to handle zeros better.  It still isn't right because the
compiler thinks that +0.0d0 and -0.0d0 are similar as a constant.
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.12 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.12 1993/05/08 04:49:36 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     (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 ram 1.6 ;;; 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 wlott 1.1 (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 ram 1.6 ;;; 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 wlott 1.1 (defun expt (base power)
133 wlott 1.3 "Returns BASE raised to the POWER."
134 wlott 1.1 (if (zerop power)
135 ram 1.6 (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 ram 1.8 (((foreach single-float double-float) rational)
155 ram 1.6 (real-expt base power '(dispatch-type base)))
156 ram 1.9 (((foreach fixnum (or bignum ratio) single-float)
157 ram 1.6 (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 wlott 1.1
175     (defun log (number &optional (base nil base-p))
176 wlott 1.3 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
177 wlott 1.1 (if base-p
178     (/ (log number) (log base))
179     (number-dispatch ((number number))
180 wlott 1.3 (((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 wlott 1.1
190     (defun sqrt (number)
191     "Return the square root of NUMBER."
192     (number-dispatch ((number number))
193 wlott 1.3 (((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 wlott 1.1
203    
204     ;;;; Trigonometic and Related Functions
205    
206 wlott 1.2 (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 wlott 1.1
224     (defun phase (number)
225     "Returns the angle part of the polar representation of a complex number.
226 wlott 1.3 For complex numbers, this is (atan (imagpart number) (realpart number)).
227 wlott 1.1 For non-complex positive numbers, this is 0. For non-complex negative
228     numbers this is PI."
229 wlott 1.3 (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 wlott 1.1
241    
242     (defun sin (number)
243 wlott 1.3 "Return the sine of NUMBER."
244 wlott 1.1 (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 wlott 1.3 "Return the cosine of NUMBER."
253 wlott 1.1 (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 wlott 1.3 "Return the tangent of NUMBER."
262 wlott 1.1 (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 wlott 1.2 (defun cis (theta)
271 wlott 1.3 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
272 wlott 1.2 (if (complexp theta)
273     (error "Argument to CIS is complex: ~S" theta)
274     (complex (cos theta) (sin theta))))
275 wlott 1.1
276 wlott 1.3 (proclaim '(inline mult-by-i))
277     (defun mult-by-i (number)
278     (complex (imagpart number)
279     (- (realpart number))))
280 wlott 1.1
281 wlott 1.3 (defun complex-asin (number)
282     (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
283    
284 wlott 1.1 (defun asin (number)
285 wlott 1.3 "Return the arc sine of NUMBER."
286 wlott 1.1 (number-dispatch ((number number))
287 wlott 1.3 ((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 wlott 1.1 ((complex)
298 wlott 1.3 (complex-asin number))))
299 wlott 1.1
300 wlott 1.3 (defun complex-acos (number)
301     (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))
302 wlott 1.1
303     (defun acos (number)
304 wlott 1.3 "Return the arc cosine of NUMBER."
305 wlott 1.1 (number-dispatch ((number number))
306 wlott 1.3 ((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 wlott 1.1 ((complex)
317 wlott 1.3 (complex-acos number))))
318 wlott 1.1
319    
320     (defun atan (y &optional (x nil xp))
321 wlott 1.3 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
322 wlott 1.1 (if xp
323 wlott 1.12 (flet ((atan2 (y x)
324     (declare (type double-float y x))
325     (if (plusp (float-sign x))
326     ;; X is either +0.0 or > 0.0
327     (if (zerop x)
328     ;; X is +0.0.
329     (if (plusp (float-sign y))
330     ;; Y is either +0.0 or > 0.0
331     (if (zerop y) +0.0d0 (/ pi 2))
332     ;; Y is either -0.0 or < 0.0
333     (if (zerop y) -0.0d0 (/ pi 2)))
334     ;; X is > 0.0
335     (if (plusp (float-sign y))
336     ;; Y is either +0.0 or > 0.0
337     (if (zerop y) +0.0d0 (%atan2 y x))
338     ;; Y is either -0.0 or < 0.0
339     (if (zerop y) -0.0d0 (%atan2 y x))))
340     ;; X is either -0.0 or < 0.0
341     (if (zerop x)
342     ;; X is -0.0
343     (if (plusp (float-sign y))
344     ;; Y is either +0.0 or > 0.0
345     (if (zerop y) pi (/ pi 2))
346     ;; Y is either -0.0 or < 0.0
347     (if (zerop y) (- pi) (- (/ pi 2))))
348     ;; X is < 0.0
349     (if (zerop y) pi (%atan2 y x))))))
350     (atan2 (number-dispatch ((y real))
351     (((foreach fixnum bignum ratio single-float))
352     (coerce y 'double-float))
353     ((double-float) y))
354     (number-dispatch ((x real))
355     (((foreach fixnum bignum ratio single-float))
356     (coerce x 'double-float))
357     ((double-float) x))))
358 wlott 1.1 (number-dispatch ((y number))
359     (handle-reals %atan y)
360     ((complex)
361 wlott 1.3 (let ((im (imagpart y))
362     (re (realpart y)))
363     (/ (- (log (complex (- 1 im) re))
364     (log (complex (+ 1 im) (- re))))
365     (complex 0 2)))))))
366    
367    
368     (defun sinh (number)
369     "Return the hyperbolic sine of NUMBER."
370     (/ (- (exp number) (exp (- number))) 2))
371    
372     (defun cosh (number)
373     "Return the hyperbolic cosine of NUMBER."
374     (/ (+ (exp number) (exp (- number))) 2))
375    
376     (defun tanh (number)
377     "Return the hyperbolic tangent of NUMBER."
378     (/ (- (exp number) (exp (- number)))
379     (+ (exp number) (exp (- number)))))
380    
381    
382     (defun asinh (number)
383     "Return the hyperbolic arc sine of NUMBER."
384     (log (+ number (sqrt (1+ (* number number))))))
385    
386     (defun acosh (number)
387     "Return the hyperbolic arc cosine of NUMBER."
388     (log (+ number (* (1+ number) (sqrt (/ (1- number) (1+ number)))))))
389    
390     (defun atanh (number)
391     "Return the hyperbolic arc tangent of NUMBER."
392     (log (* (1+ number) (sqrt (/ (- 1 (* number number)))))))
393    

  ViewVC Help
Powered by ViewVC 1.1.5