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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5