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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations)
Fri Jul 20 23:52:13 1990 UTC (23 years, 9 months ago) by wlott
Branch: MAIN
Initial revision
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 ;;; $Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.1 1990/07/20 23:52:13 wlott Exp $
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 (defconstant e 2.71828182845904523536028747135266249775724709369996L0)
25
26 (defmacro def-math-rtn (name num-args)
27 (let ((function (intern (concatenate 'simple-string
28 "%"
29 (string-upcase name)))))
30 `(def-c-routine (,name ,function) (double-float)
31 ,@(let ((results nil))
32 (dotimes (i num-args (nreverse results))
33 (push (list (intern (format nil "ARG-~D" i))
34 'double-float)
35 results))))))
36
37 (eval-when (compile load eval)
38
39 (defun handle-reals (function var)
40 `((((foreach fixnum single-float bignum ratio))
41 (coerce (,function (coerce ,var 'double-float)) 'single-float))
42 ((double-float)
43 (,function ,var))))
44
45 ); eval-when (compile load eval)
46
47
48 ;;;; Stubs for the Unix math library.
49
50 ;;; Please refer to the Unix man pages for details about these routines.
51
52 ;;; Trigonometric.
53 (def-math-rtn "sin" 1)
54 (def-math-rtn "cos" 1)
55 (def-math-rtn "tan" 1)
56 (def-math-rtn "asin" 1)
57 (def-math-rtn "acos" 1)
58 (def-math-rtn "atan" 1)
59 (def-math-rtn "atan2" 2)
60 (def-math-rtn "sinh" 1)
61 (def-math-rtn "cosh" 1)
62 (def-math-rtn "tanh" 1)
63 (def-math-rtn "asinh" 1)
64 (def-math-rtn "acosh" 1)
65 (def-math-rtn "atanh" 1)
66
67 ;;; Exponential and Logarithmic.
68 (def-math-rtn "exp" 1)
69 (def-math-rtn "expm1" 1)
70 (def-math-rtn "log" 1)
71 (def-math-rtn "log10" 1)
72 (def-math-rtn "log1p" 1)
73 (def-math-rtn "pow" 2)
74 (def-math-rtn "cbrt" 1)
75 (def-math-rtn "sqrt" 1)
76 (def-math-rtn "hypot" 2)
77
78
79 ;;;; Power functions.
80
81 (defun exp (number)
82 "Return e raised to the power NUMBER."
83 (number-dispatch ((number number))
84 (handle-reals %exp number)
85 ((complex)
86 (* (exp (realpart number))
87 (cis (imagpart number))))))
88
89 ;;; INTEXP -- Handle the rational base, integer power case.
90
91 (defparameter *intexp-maximum-exponent* 10000)
92
93 (defun intexp (base power)
94 (when (> (abs power) *intexp-maximum-exponent*)
95 (cerror "Continue with calculation."
96 "The absolute value of ~S exceeds ~S."
97 power '*intexp-maximum-exponent* base power))
98 (cond ((minusp power)
99 (/ (intexp base (- power))))
100 ((eql base 2)
101 (ash 1 power))
102 (t
103 (do ((nextn (ash power -1) (ash power -1))
104 (total (if (oddp power) base 1)
105 (if (oddp power) (* base total) total)))
106 ((zerop nextn) total)
107 (setq base (* base base))
108 (setq power nextn)))))
109
110 ;;; This function calculates x raised to the nth power. It separates
111 ;;; the cases by the type of n, for efficiency reasons, as powers can
112 ;;; be calculated more efficiently if n is a positive integer, Therefore,
113 ;;; All integer values of n are calculated as positive integers, and
114 ;;; inverted if negative.
115
116 (defun expt (base power)
117 "Returns x raised to the nth power."
118 (if (zerop power)
119 ;; This is wrong if power isn't an integer.
120 (typecase (realpart base)
121 (single-float (coerce 1 'single-float))
122 (double-float (coerce 1 'double-float))
123 (t 1))
124 (number-dispatch ((base number) (power number))
125 (((foreach fixnum bignum ratio (complex rational)) integer)
126 (intexp base power))
127 (((foreach single-float double-float) integer)
128 (coerce (%pow (coerce base 'double-float)
129 (coerce power 'double-float))
130 '(dispatch-type base)))
131 (((foreach fixnum bignum ratio single-float)
132 (foreach ratio single-float))
133 (coerce (%pow (coerce base 'double-float)
134 (coerce power 'double-float))
135 'single-float))
136 (((foreach fixnum bignum ratio single-float double-float) double-float)
137 (%pow (coerce base 'double-float) (coerce power 'double-float)))
138 (((complex rational) ratio)
139 (* (expt (abs base) power)
140 (cis (* power (phase base)))))
141 (((complex float) (foreach integer ratio))
142 (* (expt (abs base) power)
143 (cis (* power (phase base)))))
144 (((foreach fixnum bignum ratio single-float double-float) complex)
145 (if (minusp base)
146 (/ (exp (* power (log (- base)))))
147 (exp (* power (log base)))))
148 (((foreach (complex float) (complex rational)) complex)
149 (exp (* power (log base)))))))
150
151 (defun log (number &optional (base nil base-p))
152 (if base-p
153 (/ (log number) (log base))
154 (number-dispatch ((number number))
155 (handle-reals %log number)
156 ((complex)
157 (complex (log (abs number))
158 (phase number))))))
159
160 (defun sqrt (number)
161 "Return the square root of NUMBER."
162 (number-dispatch ((number number))
163 (handle-reals %sqrt number)
164 ((complex)
165 (* (exp (/ (complex 0 (phase number)) 2))
166 (sqrt (abs number))))))
167
168 ;;; ISQRT: Integer square root - isqrt(n)**2 <= n
169 ;;; Upper and lower bounds on the result are estimated using integer-length.
170 ;;; On each iteration, one of the bounds is replaced by their mean.
171 ;;; The lower bound is returned when the bounds meet or differ by only 1.
172 ;;; Initial bounds guarantee that lg(sqrt(n)) = lg(n)/2 iterations suffice.
173
174 (defun isqrt (n)
175 "Returns the root of the nearest integer less than
176 n which is a perfect square."
177 (if (and (integerp n) (not (minusp n)))
178 (do* ((lg (integer-length n))
179 (lo (ash 1 (ash (1- lg) -1)))
180 (hi (+ lo (ash lo (if (oddp lg) -1 0))))) ;tighten by 3/4 if possible.
181 ((<= (1- hi) lo) lo)
182 (let ((mid (ash (+ lo hi) -1)))
183 (if (<= (* mid mid) n) (setq lo mid) (setq hi mid))))
184 (error "Isqrt: ~S argument must be a nonnegative integer" n)))
185
186
187
188 ;;;; Trigonometic and Related Functions
189
190 ;;; ABS is in numbers.lisp
191
192 (defun phase (number)
193 "Returns the angle part of the polar representation of a complex number.
194 For non-complex positive numbers, this is 0. For non-complex negative
195 numbers this is PI."
196 (if (zerop number)
197 (if (typep number 'double-float) 0.0d0 0.0)
198 (atan (imagpart number) (realpart number))))
199
200 #|
201 (if (complexp number)
202 (let ((ipart (imagpart number))
203 (rpart (realpart number)))
204 (if (zerop rpart)
205 (if (minusp ipart)
206 (if (long-float-p ipart) (- %long-pi/2) (- %short-pi/2))
207 (if (long-float-p ipart) %long-pi/2 %short-pi/2))
208 (atan (/ (imagpart number) (realpart number)))))
209 (if (minusp number)
210 (if (long-float-p number) pi %short-pi)
211 (if (long-float-p number) 0.0l0 0.0))))
212 |#
213
214 (defun sin (number)
215 (number-dispatch ((number number))
216 (handle-reals %sin number)
217 ((complex)
218 (let ((x (realpart number))
219 (y (imagpart number)))
220 (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))))
221
222 (defun cos (number)
223 (number-dispatch ((number number))
224 (handle-reals %cos number)
225 ((complex)
226 (let ((x (realpart number))
227 (y (imagpart number)))
228 (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))))
229
230 (defun tan (number)
231 (number-dispatch ((number number))
232 (handle-reals %tan number)
233 ((complex)
234 (let* ((num (sin number))
235 (denom (cos number)))
236 (if (zerop denom) (error "~S undefined tangent." number)
237 (/ num denom))))))
238
239 ;;; CIS is in numbers.lisp
240
241 #+nil
242 (defun in-asin-domain (z)
243 (or (< (- (/ pi 2.0d0)) (realpart z) (/ pi 2.0d0))
244 (and (= (realpart z) (- (/ pi 2.0d0)))
245 (>= (imagpart z) 0))
246 (and (= (realpart z) (/ pi 2.0d0))
247 (<= (imagpart z) 0))))
248
249 (defun asin (number)
250 (number-dispatch ((number number))
251 (handle-reals %asin number)
252 ((complex)
253 (error "Can't hack complex ASIN yet: ~S" number))))
254
255 #+nil
256 (defun in-acos-domain (z)
257 (or (< 0 (realpart z) pi)
258 (and (= 0 (realpart z))
259 (>= (imagpart z) 0))
260 (and (= (realpart z) pi)
261 (<= (imagpart z) 0))))
262
263 (defun acos (number)
264 (number-dispatch ((number number))
265 (handle-reals %acos number)
266 ((complex)
267 (error "Can't hack complex ACOS yet: ~S" number))))
268
269
270 (defun atan (y &optional (x nil xp))
271 (if xp
272 (if (and (zerop x) (zerop y))
273 (error "Both args to ATAN can't be zero.")
274 (number-dispatch ((y real) (x real))
275 (((foreach fixnum bignum ratio single-float)
276 (foreach fixnum bignum ratio single-float))
277 (coerce (%atan2 (coerce y 'double-float)
278 (coerce x 'double-float))
279 'single-float))
280 ((double-float (foreach fixnum bignum ratio single-float))
281 (%atan2 y (coerce x 'double-float)))
282 (((foreach fixnum bignum ratio single-float double-float)
283 double-float)
284 (%atan2 (coerce y 'double-float) x))))
285 (number-dispatch ((y number))
286 (handle-reals %atan y)
287 ((complex)
288 (error "Can't handle complex ATAN yet: ~S" y)))))

  ViewVC Help
Powered by ViewVC 1.1.5