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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.9.1.1 - (hide annotations) (vendor branch)
Fri Jan 24 08:47:04 1992 UTC (22 years, 2 months ago) by wlott
Changes since 1.9: +2 -2 lines
Mods for new aliens.
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.9.1.1 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.9.1.1 1992/01/24 08:47:04 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.9.1.1 (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     ;;; ISQRT: Integer square root - isqrt(n)**2 <= n
204     ;;; Upper and lower bounds on the result are estimated using integer-length.
205     ;;; On each iteration, one of the bounds is replaced by their mean.
206     ;;; The lower bound is returned when the bounds meet or differ by only 1.
207     ;;; Initial bounds guarantee that lg(sqrt(n)) = lg(n)/2 iterations suffice.
208    
209     (defun isqrt (n)
210     "Returns the root of the nearest integer less than
211     n which is a perfect square."
212     (if (and (integerp n) (not (minusp n)))
213     (do* ((lg (integer-length n))
214     (lo (ash 1 (ash (1- lg) -1)))
215     (hi (+ lo (ash lo (if (oddp lg) -1 0))))) ;tighten by 3/4 if possible.
216     ((<= (1- hi) lo) lo)
217     (let ((mid (ash (+ lo hi) -1)))
218     (if (<= (* mid mid) n) (setq lo mid) (setq hi mid))))
219     (error "Isqrt: ~S argument must be a nonnegative integer" n)))
220    
221    
222    
223     ;;;; Trigonometic and Related Functions
224    
225 wlott 1.2 (defun abs (number)
226     "Returns the absolute value of the number."
227     (number-dispatch ((number number))
228     (((foreach single-float double-float fixnum rational))
229     (abs number))
230     ((complex)
231     (let ((rx (realpart number))
232     (ix (imagpart number)))
233     (etypecase rx
234     (rational
235     (sqrt (+ (* rx rx) (* ix ix))))
236     (single-float
237     (coerce (%hypot (coerce rx 'double-float)
238     (coerce ix 'double-float))
239     'single-float))
240     (double-float
241     (%hypot rx ix)))))))
242 wlott 1.1
243     (defun phase (number)
244     "Returns the angle part of the polar representation of a complex number.
245 wlott 1.3 For complex numbers, this is (atan (imagpart number) (realpart number)).
246 wlott 1.1 For non-complex positive numbers, this is 0. For non-complex negative
247     numbers this is PI."
248 wlott 1.3 (etypecase number
249     ((or rational single-float)
250     (if (minusp number)
251     (coerce pi 'single-float)
252     0.0f0))
253     (double-float
254     (if (minusp number)
255     (coerce pi 'double-float)
256     0.0d0))
257     (complex
258     (atan (imagpart number) (realpart number)))))
259 wlott 1.1
260    
261     (defun sin (number)
262 wlott 1.3 "Return the sine of NUMBER."
263 wlott 1.1 (number-dispatch ((number number))
264     (handle-reals %sin number)
265     ((complex)
266     (let ((x (realpart number))
267     (y (imagpart number)))
268     (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))))
269    
270     (defun cos (number)
271 wlott 1.3 "Return the cosine of NUMBER."
272 wlott 1.1 (number-dispatch ((number number))
273     (handle-reals %cos number)
274     ((complex)
275     (let ((x (realpart number))
276     (y (imagpart number)))
277     (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))))
278    
279     (defun tan (number)
280 wlott 1.3 "Return the tangent of NUMBER."
281 wlott 1.1 (number-dispatch ((number number))
282     (handle-reals %tan number)
283     ((complex)
284     (let* ((num (sin number))
285     (denom (cos number)))
286     (if (zerop denom) (error "~S undefined tangent." number)
287     (/ num denom))))))
288    
289 wlott 1.2 (defun cis (theta)
290 wlott 1.3 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
291 wlott 1.2 (if (complexp theta)
292     (error "Argument to CIS is complex: ~S" theta)
293     (complex (cos theta) (sin theta))))
294 wlott 1.1
295 wlott 1.3 (proclaim '(inline mult-by-i))
296     (defun mult-by-i (number)
297     (complex (imagpart number)
298     (- (realpart number))))
299 wlott 1.1
300 wlott 1.3 (defun complex-asin (number)
301     (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
302    
303 wlott 1.1 (defun asin (number)
304 wlott 1.3 "Return the arc sine of NUMBER."
305 wlott 1.1 (number-dispatch ((number number))
306 wlott 1.3 ((rational)
307     (if (or (> number 1) (< number -1))
308     (complex-asin number)
309     (coerce (%asin (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-asin number)
314     (coerce (%asin (coerce number 'double-float))
315     '(dispatch-type number))))
316 wlott 1.1 ((complex)
317 wlott 1.3 (complex-asin number))))
318 wlott 1.1
319 wlott 1.3 (defun complex-acos (number)
320     (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))
321 wlott 1.1
322     (defun acos (number)
323 wlott 1.3 "Return the arc cosine of NUMBER."
324 wlott 1.1 (number-dispatch ((number number))
325 wlott 1.3 ((rational)
326     (if (or (> number 1) (< number -1))
327     (complex-acos number)
328     (coerce (%acos (coerce number 'double-float)) 'single-float)))
329     (((foreach single-float double-float))
330     (if (or (> number (coerce 1 '(dispatch-type number)))
331     (< number (coerce -1 '(dispatch-type number))))
332     (complex-acos number)
333     (coerce (%acos (coerce number 'double-float))
334     '(dispatch-type number))))
335 wlott 1.1 ((complex)
336 wlott 1.3 (complex-acos number))))
337 wlott 1.1
338    
339     (defun atan (y &optional (x nil xp))
340 wlott 1.3 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
341 wlott 1.1 (if xp
342     (if (and (zerop x) (zerop y))
343 ram 1.4 (if (plusp (float-sign (float x)))
344     y
345     (if (minusp (float-sign (float y)))
346     (- pi)
347     pi))
348 wlott 1.1 (number-dispatch ((y real) (x real))
349     (((foreach fixnum bignum ratio single-float)
350     (foreach fixnum bignum ratio single-float))
351     (coerce (%atan2 (coerce y 'double-float)
352     (coerce x 'double-float))
353     'single-float))
354     ((double-float (foreach fixnum bignum ratio single-float))
355     (%atan2 y (coerce x 'double-float)))
356     (((foreach fixnum bignum ratio single-float double-float)
357     double-float)
358     (%atan2 (coerce y 'double-float) x))))
359     (number-dispatch ((y number))
360     (handle-reals %atan y)
361     ((complex)
362 wlott 1.3 (let ((im (imagpart y))
363     (re (realpart y)))
364     (/ (- (log (complex (- 1 im) re))
365     (log (complex (+ 1 im) (- re))))
366     (complex 0 2)))))))
367    
368    
369     (defun sinh (number)
370     "Return the hyperbolic sine of NUMBER."
371     (/ (- (exp number) (exp (- number))) 2))
372    
373     (defun cosh (number)
374     "Return the hyperbolic cosine of NUMBER."
375     (/ (+ (exp number) (exp (- number))) 2))
376    
377     (defun tanh (number)
378     "Return the hyperbolic tangent of NUMBER."
379     (/ (- (exp number) (exp (- number)))
380     (+ (exp number) (exp (- number)))))
381    
382    
383     (defun asinh (number)
384     "Return the hyperbolic arc sine of NUMBER."
385     (log (+ number (sqrt (1+ (* number number))))))
386    
387     (defun acosh (number)
388     "Return the hyperbolic arc cosine of NUMBER."
389     (log (+ number (* (1+ number) (sqrt (/ (1- number) (1+ number)))))))
390    
391     (defun atanh (number)
392     "Return the hyperbolic arc tangent of NUMBER."
393     (log (* (1+ number) (sqrt (/ (- 1 (* number number)))))))
394    

  ViewVC Help
Powered by ViewVC 1.1.5