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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.26 - (show annotations)
Thu Feb 19 03:49:49 1998 UTC (16 years, 2 months ago) by dtc
Branch: MAIN
Changes since 1.25: +4 -4 lines
Fix float-trapping-nan-p which was returning T for quiet NaN and Nil
of trapping NaN.
1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2 ;;;
3 ;;; **********************************************************************
4 ;;; 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 ;;;
7 (ext:file-comment
8 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.26 1998/02/19 03:49:49 dtc Exp $")
9 ;;;
10 ;;; **********************************************************************
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 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
27 ;;; call, and saves number consing to boot.
28 ;;;
29 (defmacro def-math-rtn (name num-args)
30 (let ((function (intern (concatenate 'simple-string
31 "%"
32 (string-upcase name)))))
33 `(progn
34 (proclaim '(inline ,function))
35 (export ',function)
36 (alien:def-alien-routine (,name ,function) double-float
37 ,@(let ((results nil))
38 (dotimes (i num-args (nreverse results))
39 (push (list (intern (format nil "ARG-~D" i))
40 'double-float)
41 results)))))))
42
43 (eval-when (compile load eval)
44
45 (defun handle-reals (function var)
46 `((((foreach fixnum single-float bignum ratio))
47 (coerce (,function (coerce ,var 'double-float)) 'single-float))
48 ((double-float)
49 (,function ,var))))
50
51 ); eval-when (compile load eval)
52
53
54 ;;;; Stubs for the Unix math library.
55
56 ;;; Please refer to the Unix man pages for details about these routines.
57
58 ;;; Trigonometric.
59 #-x86 (def-math-rtn "sin" 1)
60 #-x86 (def-math-rtn "cos" 1)
61 #-x86 (def-math-rtn "tan" 1)
62 (def-math-rtn "asin" 1)
63 (def-math-rtn "acos" 1)
64 #-x86 (def-math-rtn "atan" 1)
65 #-x86 (def-math-rtn "atan2" 2)
66 (def-math-rtn "sinh" 1)
67 (def-math-rtn "cosh" 1)
68 (def-math-rtn "tanh" 1)
69 (def-math-rtn "asinh" 1)
70 (def-math-rtn "acosh" 1)
71 (def-math-rtn "atanh" 1)
72
73 ;;; Exponential and Logarithmic.
74 #-x86 (def-math-rtn "exp" 1)
75 #-x86 (def-math-rtn "log" 1)
76 #-x86 (def-math-rtn "log10" 1)
77 (def-math-rtn "pow" 2)
78 #-x86 (def-math-rtn "sqrt" 1)
79 (def-math-rtn "hypot" 2)
80 #-hpux
81 (def-math-rtn "log1p" 1)
82
83 #+x86 ;; These are needed for use by byte-compiled files.
84 (progn
85 (defun %sin (x)
86 (declare (double-float x)
87 (values double-float))
88 (%sin x))
89 (defun %sin-quick (x)
90 (declare (double-float x)
91 (values double-float))
92 (%sin-quick x))
93 (defun %cos (x)
94 (declare (double-float x)
95 (values double-float))
96 (%cos x))
97 (defun %cos-quick (x)
98 (declare (double-float x)
99 (values double-float))
100 (%cos-quick x))
101 (defun %tan (x)
102 (declare (double-float x)
103 (values double-float))
104 (%tan x))
105 (defun %tan-quick (x)
106 (declare (double-float x)
107 (values double-float))
108 (%tan-quick x))
109 (defun %atan (x)
110 (declare (double-float x)
111 (values double-float))
112 (%atan x))
113 (defun %atan2 (x y)
114 (declare (double-float x y)
115 (values double-float))
116 (%atan2 x y))
117 (defun %exp (x)
118 (declare (double-float x)
119 (values double-float))
120 (%exp x))
121 (defun %log (x)
122 (declare (double-float x)
123 (values double-float))
124 (%log x))
125 (defun %log10 (x)
126 (declare (double-float x)
127 (values double-float))
128 (%log10 x))
129 #+nil ;; notyet
130 (defun %pow (x y)
131 (declare (type (double-float 0d0) x)
132 (double-float y)
133 (values (double-float 0d0)))
134 (%pow x y))
135 (defun %sqrt (x)
136 (declare (double-float x)
137 (values double-float))
138 (%sqrt x))
139 (defun %scalbn (f ex)
140 (declare (double-float f)
141 (type (signed-byte 32) ex)
142 (values double-float))
143 (%scalbn f ex))
144 (defun %scalb (f ex)
145 (declare (double-float f ex)
146 (values double-float))
147 (%scalb f ex))
148 (defun %logb (x)
149 (declare (double-float x)
150 (values double-float))
151 (%logb x))
152 ) ; progn
153
154
155 ;;;; Power functions.
156
157 (defun exp (number)
158 "Return e raised to the power NUMBER."
159 (number-dispatch ((number number))
160 (handle-reals %exp number)
161 ((complex)
162 (* (exp (realpart number))
163 (cis (imagpart number))))))
164
165 ;;; INTEXP -- Handle the rational base, integer power case.
166
167 (defparameter *intexp-maximum-exponent* 10000)
168
169 ;;; This function precisely calculates base raised to an integral power. It
170 ;;; separates the cases by the sign of power, for efficiency reasons, as powers
171 ;;; can be calculated more efficiently if power is a positive integer. Values
172 ;;; of power are calculated as positive integers, and inverted if negative.
173 ;;;
174 (defun intexp (base power)
175 (when (> (abs power) *intexp-maximum-exponent*)
176 (cerror "Continue with calculation."
177 "The absolute value of ~S exceeds ~S."
178 power '*intexp-maximum-exponent* base power))
179 (cond ((minusp power)
180 (/ (intexp base (- power))))
181 ((eql base 2)
182 (ash 1 power))
183 (t
184 (do ((nextn (ash power -1) (ash power -1))
185 (total (if (oddp power) base 1)
186 (if (oddp power) (* base total) total)))
187 ((zerop nextn) total)
188 (setq base (* base base))
189 (setq power nextn)))))
190
191
192 ;;; EXPT -- Public
193 ;;;
194 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
195 ;;; floating point stuff. If both args are real, we try %POW right off,
196 ;;; assuming it will return 0 if the result may be complex. If so, we call
197 ;;; COMPLEX-POW which directly computes the complex result. We also separate
198 ;;; the complex-real and real-complex cases from the general complex case.
199 ;;;
200 (defun expt (base power)
201 "Returns BASE raised to the POWER."
202 (if (zerop power)
203 (1+ (* base power))
204 (labels (;; determine if the double float is an integer.
205 ;; 0 - not an integer
206 ;; 1 - an odd int
207 ;; 2 - an even int
208 (isint (ihi lo)
209 (declare (type (unsigned-byte 31) ihi)
210 (type (unsigned-byte 32) lo)
211 (optimize (speed 3) (safety 0)))
212 (let ((isint 0))
213 (declare (type fixnum isint))
214 (cond ((>= ihi #x43400000) ; exponent >= 53
215 (setq isint 2))
216 ((>= ihi #x3ff00000)
217 (let ((k (- (ash ihi -20) #x3ff))) ; exponent
218 (declare (type (mod 53) k))
219 (cond ((> k 20)
220 (let* ((shift (- 52 k))
221 (j (logand (ash lo (- shift))))
222 (j2 (ash j shift)))
223 (declare (type (mod 32) shift)
224 (type (unsigned-byte 32) j j2))
225 (when (= j2 lo)
226 (setq isint (- 2 (logand j 1))))))
227 ((= lo 0)
228 (let* ((shift (- 20 k))
229 (j (ash ihi (- shift)))
230 (j2 (ash j shift)))
231 (declare (type (mod 32) shift)
232 (type (unsigned-byte 31) j j2))
233 (when (= j2 ihi)
234 (setq isint (- 2 (logand j 1))))))))))
235 isint))
236 (real-expt (x y rtype)
237 (let ((x (coerce x 'double-float))
238 (y (coerce y 'double-float)))
239 (declare (double-float x y))
240 (let* ((x-hi (kernel:double-float-high-bits x))
241 (x-lo (kernel:double-float-low-bits x))
242 (x-ihi (logand x-hi #x7fffffff))
243 (y-hi (kernel:double-float-high-bits y))
244 (y-lo (kernel:double-float-low-bits y))
245 (y-ihi (logand y-hi #x7fffffff)))
246 (declare (type (signed-byte 32) x-hi y-hi)
247 (type (unsigned-byte 31) x-ihi y-ihi)
248 (type (unsigned-byte 32) x-lo y-lo))
249 ;; y==zero: x**0 = 1
250 (when (zerop (logior y-ihi y-lo))
251 (return-from real-expt (coerce 1d0 rtype)))
252 ;; +-NaN return x+y
253 (when (or (> x-ihi #x7ff00000)
254 (and (= x-ihi #x7ff00000) (/= x-lo 0))
255 (> y-ihi #x7ff00000)
256 (and (= y-ihi #x7ff00000) (/= y-lo 0)))
257 (return-from real-expt (coerce (+ x y) rtype)))
258 (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
259 (declare (type fixnum yisint))
260 ;; special value of y
261 (when (and (zerop y-lo) (= y-ihi #x7ff00000))
262 ;; y is +-inf
263 (return-from real-expt
264 (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
265 ;; +-1**inf is NaN
266 (coerce (- y y) rtype))
267 ((>= x-ihi #x3ff00000)
268 ;; (|x|>1)**+-inf = inf,0
269 (if (>= y-hi 0)
270 (coerce y rtype)
271 (coerce 0 rtype)))
272 (t
273 ;; (|x|<1)**-,+inf = inf,0
274 (if (< y-hi 0)
275 (coerce (- y) rtype)
276 (coerce 0 rtype))))))
277
278 (let ((abs-x (abs x)))
279 (declare (double-float abs-x))
280 ;; special value of x
281 (when (and (zerop x-lo)
282 (or (= x-ihi #x7ff00000) (zerop x-ihi)
283 (= x-ihi #x3ff00000)))
284 ;; x is +-0,+-inf,+-1
285 (let ((z (if (< y-hi 0)
286 (/ 1 abs-x) ; z = (1/|x|)
287 abs-x)))
288 (declare (double-float z))
289 (when (< x-hi 0)
290 (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
291 ;; (-1)**non-int
292 (let ((y*pi (* y pi)))
293 (declare (double-float y*pi))
294 (return-from real-expt
295 (complex
296 (coerce (%cos y*pi) rtype)
297 (coerce (%sin y*pi) rtype)))))
298 ((= yisint 1)
299 ;; (x<0)**odd = -(|x|**odd)
300 (setq z (- z)))))
301 (return-from real-expt (coerce z rtype))))
302
303 (if (>= x-hi 0)
304 ;; x>0
305 (coerce (kernel::%pow x y) rtype)
306 ;; x<0
307 (let ((pow (kernel::%pow abs-x y)))
308 (declare (double-float pow))
309 (case yisint
310 (1 ; Odd
311 (coerce (* -1d0 pow) rtype))
312 (2 ; Even
313 (coerce pow rtype))
314 (t ; Non-integer
315 (let ((y*pi (* y pi)))
316 (declare (double-float y*pi))
317 (complex
318 (coerce (* pow (%cos y*pi)) rtype)
319 (coerce (* pow (%sin y*pi)) rtype)))))))))))))
320 (declare (inline real-expt))
321 (number-dispatch ((base number) (power number))
322 (((foreach fixnum (or bignum ratio) (complex rational)) integer)
323 (intexp base power))
324 (((foreach single-float double-float) rational)
325 (real-expt base power '(dispatch-type base)))
326 (((foreach fixnum (or bignum ratio) single-float)
327 (foreach ratio single-float))
328 (real-expt base power 'single-float))
329 (((foreach fixnum (or bignum ratio) single-float double-float)
330 double-float)
331 (real-expt base power 'double-float))
332 ((double-float single-float)
333 (real-expt base power 'double-float))
334 (((foreach (complex rational) (complex float)) rational)
335 (* (expt (abs base) power)
336 (cis (* power (phase base)))))
337 (((foreach fixnum (or bignum ratio) single-float double-float)
338 complex)
339 (if (and (zerop base) (plusp (realpart power)))
340 (* base power)
341 (exp (* power (log base)))))
342 (((foreach (complex float) (complex rational))
343 (foreach complex double-float single-float))
344 (if (and (zerop base) (plusp (realpart power)))
345 (* base power)
346 (exp (* power (log base)))))))))
347
348 (defun log (number &optional (base nil base-p))
349 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
350 (if base-p
351 (if (zerop base)
352 base ; ANSI spec
353 (/ (log number) (log base)))
354 (number-dispatch ((number number))
355 (((foreach fixnum bignum ratio))
356 (if (minusp number)
357 (complex (log (- number)) (coerce pi 'single-float))
358 (coerce (%log (coerce number 'double-float)) 'single-float)))
359 (((foreach single-float double-float))
360 ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
361 ;; Since this doesn't seem to be an implementation issue
362 ;; I (pw) take the Kahan result.
363 (if (< (float-sign number)
364 (coerce 0 '(dispatch-type number)))
365 (complex (log (- number)) (coerce pi '(dispatch-type number)))
366 (coerce (%log (coerce number 'double-float))
367 '(dispatch-type number))))
368 ((complex)
369 (complex-log number)))))
370
371 (defun sqrt (number)
372 "Return the square root of NUMBER."
373 (number-dispatch ((number number))
374 (((foreach fixnum bignum ratio))
375 (if (minusp number)
376 (complex-sqrt number)
377 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
378 (((foreach single-float double-float))
379 (if (< (float-sign number)
380 (coerce 0 '(dispatch-type number)))
381 (complex-sqrt number)
382 (coerce (%sqrt (coerce number 'double-float))
383 '(dispatch-type number))))
384 ((complex)
385 (complex-sqrt number))))
386
387
388 ;;;; Trigonometic and Related Functions
389
390 (defun abs (number)
391 "Returns the absolute value of the number."
392 (number-dispatch ((number number))
393 (((foreach single-float double-float fixnum rational))
394 (abs number))
395 ((complex)
396 (let ((rx (realpart number))
397 (ix (imagpart number)))
398 (etypecase rx
399 (rational
400 (sqrt (+ (* rx rx) (* ix ix))))
401 (single-float
402 (coerce (%hypot (coerce rx 'double-float)
403 (coerce ix 'double-float))
404 'single-float))
405 (double-float
406 (%hypot rx ix)))))))
407
408 (defun phase (number)
409 "Returns the angle part of the polar representation of a complex number.
410 For complex numbers, this is (atan (imagpart number) (realpart number)).
411 For non-complex positive numbers, this is 0. For non-complex negative
412 numbers this is PI."
413 (etypecase number
414 (rational
415 (if (minusp number)
416 (coerce pi 'single-float)
417 0.0f0))
418 (single-float
419 (if (minusp (float-sign number))
420 (coerce pi 'single-float)
421 0.0f0))
422 (double-float
423 (if (minusp (float-sign number))
424 (coerce pi 'double-float)
425 0.0d0))
426 (complex
427 (atan (imagpart number) (realpart number)))))
428
429
430 (defun sin (number)
431 "Return the sine of NUMBER."
432 (number-dispatch ((number number))
433 (handle-reals %sin number)
434 ((complex)
435 (let ((x (realpart number))
436 (y (imagpart number)))
437 (complex (* (sin x) (cosh y))
438 (* (cos x) (sinh y)))))))
439
440 (defun cos (number)
441 "Return the cosine of NUMBER."
442 (number-dispatch ((number number))
443 (handle-reals %cos number)
444 ((complex)
445 (let ((x (realpart number))
446 (y (imagpart number)))
447 (complex (* (cos x) (cosh y))
448 (- (* (sin x) (sinh y))))))))
449
450 (defun tan (number)
451 "Return the tangent of NUMBER."
452 (number-dispatch ((number number))
453 (handle-reals %tan number)
454 ((complex)
455 (complex-tan number))))
456
457 (defun cis (theta)
458 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
459 (if (complexp theta)
460 (error "Argument to CIS is complex: ~S" theta)
461 (complex (cos theta) (sin theta))))
462
463 (defun asin (number)
464 "Return the arc sine of NUMBER."
465 (number-dispatch ((number number))
466 ((rational)
467 (if (or (> number 1) (< number -1))
468 (complex-asin number)
469 (coerce (%asin (coerce number 'double-float)) 'single-float)))
470 (((foreach single-float double-float))
471 (if (or (> number (coerce 1 '(dispatch-type number)))
472 (< number (coerce -1 '(dispatch-type number))))
473 (complex-asin number)
474 (coerce (%asin (coerce number 'double-float))
475 '(dispatch-type number))))
476 ((complex)
477 (complex-asin number))))
478
479 (defun acos (number)
480 "Return the arc cosine of NUMBER."
481 (number-dispatch ((number number))
482 ((rational)
483 (if (or (> number 1) (< number -1))
484 (complex-acos number)
485 (coerce (%acos (coerce number 'double-float)) 'single-float)))
486 (((foreach single-float double-float))
487 (if (or (> number (coerce 1 '(dispatch-type number)))
488 (< number (coerce -1 '(dispatch-type number))))
489 (complex-acos number)
490 (coerce (%acos (coerce number 'double-float))
491 '(dispatch-type number))))
492 ((complex)
493 (complex-acos number))))
494
495
496 (defun atan (y &optional (x nil xp))
497 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
498 (if xp
499 (flet ((atan2 (y x)
500 (declare (type double-float y x)
501 (values double-float))
502 (if (zerop x)
503 (if (zerop y)
504 (if (plusp (float-sign x))
505 y
506 (float-sign y pi))
507 (float-sign y (/ pi 2)))
508 (%atan2 y x))))
509 (number-dispatch ((y number) (x number))
510 ((double-float
511 (foreach double-float single-float fixnum bignum ratio))
512 (atan2 y (coerce x 'double-float)))
513 (((foreach single-float fixnum bignum ratio)
514 double-float)
515 (atan2 (coerce y 'double-float) x))
516 (((foreach single-float fixnum bignum ratio)
517 (foreach single-float fixnum bignum ratio))
518 (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
519 'single-float))))
520 (number-dispatch ((y number))
521 (handle-reals %atan y)
522 ((complex)
523 (complex-atan y)))))
524
525 ;; It seems that everyone has a C version of sinh, cosh, and
526 ;; tanh. Let's use these for reals because the original
527 ;; implementations based on the definitions lose big in round-off
528 ;; error. These bad definitions also mean that sin and cos for
529 ;; complex numbers can also lose big.
530
531 #+nil
532 (defun sinh (number)
533 "Return the hyperbolic sine of NUMBER."
534 (/ (- (exp number) (exp (- number))) 2))
535
536 (defun sinh (number)
537 "Return the hyperbolic sine of NUMBER."
538 (number-dispatch ((number number))
539 (handle-reals %sinh number)
540 ((complex)
541 (let ((x (realpart number))
542 (y (imagpart number)))
543 (complex (* (sinh x) (cos y))
544 (* (cosh x) (sin y)))))))
545
546 #+nil
547 (defun cosh (number)
548 "Return the hyperbolic cosine of NUMBER."
549 (/ (+ (exp number) (exp (- number))) 2))
550
551 (defun cosh (number)
552 "Return the hyperbolic cosine of NUMBER."
553 (number-dispatch ((number number))
554 (handle-reals %cosh number)
555 ((complex)
556 (let ((x (realpart number))
557 (y (imagpart number)))
558 (complex (* (cosh x) (cos y))
559 (* (sinh x) (sin y)))))))
560
561 (defun tanh (number)
562 "Return the hyperbolic tangent of NUMBER."
563 (number-dispatch ((number number))
564 (handle-reals %tanh number)
565 ((complex)
566 (complex-tanh number))))
567
568 (defun asinh (number)
569 "Return the hyperbolic arc sine of NUMBER."
570 (number-dispatch ((number number))
571 (handle-reals %asinh number)
572 ((complex)
573 (complex-asinh number))))
574
575 (defun acosh (number)
576 "Return the hyperbolic arc cosine of NUMBER."
577 (number-dispatch ((number number))
578 ((rational)
579 ;; acosh is complex if number < 1
580 (if (< number 1)
581 (complex-acosh number)
582 (coerce (%acosh (coerce number 'double-float)) 'single-float)))
583 (((foreach single-float double-float))
584 (if (< number (coerce 1 '(dispatch-type number)))
585 (complex-acosh number)
586 (coerce (%acosh (coerce number 'double-float))
587 '(dispatch-type number))))
588 ((complex)
589 (complex-acosh number))))
590
591 (defun atanh (number)
592 "Return the hyperbolic arc tangent of NUMBER."
593 (number-dispatch ((number number))
594 ((rational)
595 ;; atanh is complex if |number| > 1
596 (if (or (> number 1) (< number -1))
597 (complex-atanh number)
598 (coerce (%atanh (coerce number 'double-float)) 'single-float)))
599 (((foreach single-float double-float))
600 (if (or (> number (coerce 1 '(dispatch-type number)))
601 (< number (coerce -1 '(dispatch-type number))))
602 (complex-atanh number)
603 (coerce (%atanh (coerce number 'double-float))
604 '(dispatch-type number))))
605 ((complex)
606 (complex-atanh number))))
607
608 ;;; HP-UX does not supply a C version of log1p, so
609 ;;; use the definition.
610
611 #+hpux
612 (declaim (inline %log1p))
613 #+hpux
614 (defun %log1p (number)
615 (declare (double-float number)
616 (optimize (speed 3) (safety 0)))
617 (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
618
619
620 #+old-elfun
621 (progn
622 ;;; Here are the old definitions of the special functions, for
623 ;;; complex-valued arguments. Some of these functions suffer from
624 ;;; severe round-off error or unnecessary overflow.
625
626 (proclaim '(inline mult-by-i))
627 (defun mult-by-i (number)
628 (complex (- (imagpart number))
629 (realpart number)))
630
631 (defun complex-sqrt (x)
632 (exp (/ (log x) 2)))
633
634 (defun complex-log (x)
635 (complex (log (abs x))
636 (phase x)))
637
638 (defun complex-atanh (number)
639 (/ (- (log (1+ number)) (log (- 1 number))) 2))
640
641 (defun complex-tanh (number)
642 (/ (- (exp number) (exp (- number)))
643 (+ (exp number) (exp (- number)))))
644
645 (defun complex-acos (number)
646 (* -2 (mult-by-i (log (+ (sqrt (/ (1+ number) 2))
647 (mult-by-i (sqrt (/ (- 1 number) 2))))))))
648
649 (defun complex-acosh (number)
650 (* 2 (log (+ (sqrt (/ (1+ number) 2)) (sqrt (/ (1- number) 2))))))
651
652 (defun complex-asin (number)
653 (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
654
655 (defun complex-asinh (number)
656 (log (+ number (sqrt (1+ (* number number))))))
657
658 (defun complex-atan (y)
659 (let ((im (imagpart y))
660 (re (realpart y)))
661 (/ (- (log (complex (- 1 im) re))
662 (log (complex (+ 1 im) (- re))))
663 (complex 0 2))))
664
665 (defun complex-tan (number)
666 (let* ((num (sin number))
667 (denom (cos number)))
668 (if (zerop denom) (error "~S undefined tangent." number)
669 (/ num denom))))
670 )
671
672 #-old-specfun
673 (progn
674 ;;;;
675 ;;;; This is a set of routines that implement many elementary
676 ;;;; transcendental functions as specified by ANSI Common Lisp. The
677 ;;;; implementation is based on Kahan's paper.
678 ;;;;
679 ;;;; I believe I have accurately implemented the routines and are
680 ;;;; correct, but you may want to check for your self.
681 ;;;;
682 ;;;; These functions are written for CMU Lisp and take advantage of
683 ;;;; some of the features available there. It may be possible,
684 ;;;; however, to port this to other Lisps.
685 ;;;;
686 ;;;; Some functions are significantly more accurate than the original
687 ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
688 ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
689 ;;;; answer is pi + i*log(2-sqrt(3)).
690 ;;;;
691 ;;;; All of the implemented functions will take any number for an
692 ;;;; input, but the result will always be a either a complex
693 ;;;; single-float or a complex double-float.
694 ;;;;
695 ;;;; General functions
696 ;;;; complex-sqrt
697 ;;;; complex-log
698 ;;;; complex-atanh
699 ;;;; complex-tanh
700 ;;;; complex-acos
701 ;;;; complex-acosh
702 ;;;; complex-asin
703 ;;;; complex-asinh
704 ;;;; complex-atan
705 ;;;; complex-tan
706 ;;;;
707 ;;;; Utility functions:
708 ;;;; scalb logb
709 ;;;;
710 ;;;; Internal functions:
711 ;;;; square coerce-to-complex-type cssqs complex-log-scaled
712 ;;;;
713 ;;;;
714 ;;;; Please send any bug reports, comments, or improvements to Raymond
715 ;;;; Toy at toy@rtp.ericsson.se.
716 ;;;;
717 ;;;; References
718 ;;;;
719 ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
720 ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
721 ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
722 ;;;; Press, 1987
723 ;;;;
724 (declaim (inline square))
725 (declaim (ftype (function (double-float) (double-float 0d0)) square))
726 (defun square (x)
727 (declare (double-float x)
728 (values (double-float 0d0)))
729 (* x x))
730
731 ;; If you have these functions in libm, perhaps they should be used
732 ;; instead of these Lisp versions. These versions are probably good
733 ;; enough, especially since they are portable.
734
735 (declaim (inline scalb))
736 (defun scalb (x n)
737 "Compute 2^N * X without compute 2^N first (use properties of the
738 underlying floating-point format"
739 (declare (type double-float x)
740 (type double-float-exponent n))
741 (scale-float x n))
742
743 (defun logb (x)
744 "Compute an integer N such that 1 <= |2^N * x| < 2.
745 For the special cases, the following values are used:
746
747 x logb
748 NaN NaN
749 +/- infinity +infinity
750 0 -infinity
751 "
752 (declare (type double-float x))
753 (cond ((float-nan-p x)
754 x)
755 ((float-infinity-p x)
756 #.ext:double-float-positive-infinity)
757 ((zerop x)
758 ;; The answer is negative infinity, but we are supposed to
759 ;; signal divide-by-zero.
760 ;; (error 'division-by-zero :operation 'logb :operands (list x))
761 (/ -1.0d0 x)
762 )
763 (t
764 (multiple-value-bind (signif expon sign)
765 (decode-float x)
766 (declare (ignore signif sign))
767 ;; decode-float is almost right, except that the exponent
768 ;; is off by one
769 (1- expon)))))
770
771 ;; This function is used to create a complex number of the appropriate
772 ;; type.
773
774 (declaim (inline coerce-to-complex-type))
775 (defun coerce-to-complex-type (x y z)
776 "Create complex number with real part X and imaginary part Y such that
777 it has the same type as Z. If Z has type (complex rational), the X
778 and Y are coerced to single-float."
779 (declare (double-float x y)
780 (number z))
781 (if (subtypep (type-of (realpart z)) 'double-float)
782 (complex x y)
783 ;; Convert anything that's not a double-float to a single-float.
784 (complex (float x 1.0)
785 (float y 1.0))))
786
787 (defun cssqs (z)
788 ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
789 ;; result is r + i*k, where k is an integer.
790
791 ;; Save all FP flags
792 (let ((x (float (realpart z) 1d0))
793 (y (float (imagpart z) 1d0))
794 (k 0)
795 (rho 0d0))
796 (declare (double-float x y)
797 (type (double-float 0d0) rho)
798 (fixnum k))
799 ;; Would this be better handled using an exception handler to
800 ;; catch the overflow or underflow signal? For now, we turn all
801 ;; traps off and look at the accrued exceptions to see if any
802 ;; signal would have been raised.
803 (with-float-traps-masked (:underflow :overflow)
804 (setf rho (+ (square x) (square y)))
805 (cond ((and (or (float-nan-p rho)
806 (float-infinity-p rho))
807 (or (float-infinity-p (abs x))
808 (float-infinity-p (abs y))))
809 (setf rho #.ext:double-float-positive-infinity))
810 ((let ((threshold #.(/ least-positive-double-float
811 double-float-epsilon))
812 (traps (ldb vm::float-sticky-bits
813 (vm:floating-point-modes))))
814 ;; Overflow raised or (underflow raised and rho <
815 ;; lambda/eps)
816 (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
817 (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
818 (< rho threshold))))
819 (setf k (logb (max (abs x) (abs y))))
820 (setf rho (+ (square (scalb x (- k)))
821 (square (scalb y (- k))))))))
822 (values rho k)))
823
824 (defun complex-sqrt (z)
825 "Principle square root of Z
826
827 Z may be any number, but the result is always a complex."
828 (declare (number z))
829 (multiple-value-bind (rho k)
830 (cssqs z)
831 (declare (type (double-float 0d0) rho)
832 (fixnum k))
833 (let ((x (float (realpart z) 1.0d0))
834 (y (float (imagpart z) 1.0d0))
835 (eta 0d0)
836 (nu 0d0))
837 (declare (double-float x y eta nu))
838
839 (if (not (float-nan-p x))
840 (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
841
842 (cond ((oddp k)
843 (setf k (ash k -1)))
844 (t
845 (setf k (1- (ash k -1)))
846 (setf rho (+ rho rho))))
847
848 (setf rho (scalb (sqrt rho) k))
849
850 (setf eta rho)
851 (setf nu y)
852
853 (when (/= rho 0d0)
854 (when (not (float-infinity-p (abs nu)))
855 (setf nu (/ (/ nu rho) 2d0)))
856 (when (< x 0d0)
857 (setf eta (abs nu))
858 (setf nu (float-sign y rho))))
859 (coerce-to-complex-type eta nu z))))
860
861 (defun complex-log-scaled (z j)
862 "Compute log(2^j*z).
863
864 This is for use with J /= 0 only when |z| is huge."
865 (declare (number z)
866 (fixnum j))
867 ;; The constants t0, t1, t2 should be evaluated to machine
868 ;; precision. In addition, Kahan says the accuracy of log1p
869 ;; influences the choices of these constants but doesn't say how to
870 ;; choose them. We'll just assume his choices matches our
871 ;; implementation of log1p.
872 (let ((t0 #.(/ 1 (sqrt 2.0d0)))
873 (t1 1.2d0)
874 (t2 3d0)
875 (ln2 #.(log 2d0))
876 (x (float (realpart z) 1.0d0))
877 (y (float (imagpart z) 1.0d0)))
878 (multiple-value-bind (rho k)
879 (cssqs z)
880 (declare (type (double-float 0d0) rho)
881 (fixnum k))
882 (let ((beta (max (abs x) (abs y)))
883 (theta (min (abs x) (abs y))))
884 (declare (type (double-float 0d0) beta theta))
885 (if (and (zerop k)
886 (< t0 beta)
887 (or (<= beta t1)
888 (< rho t2)))
889 (setf rho (/ (%log1p (+ (* (- beta 1.0d0)
890 (+ beta 1.0d0))
891 (* theta theta)))
892 2d0))
893 (setf rho (+ (/ (log rho) 2d0)
894 (* (+ k j) ln2))))
895 (setf theta (atan y x))
896 (coerce-to-complex-type rho theta z)))))
897
898 (defun complex-log (z)
899 "Log of Z = log |Z| + i * arg Z
900
901 Z may be any number, but the result is always a complex."
902 (declare (number z))
903 (complex-log-scaled z 0))
904
905 ;; Let us note the following "strange" behavior. atanh 1.0d0 is
906 ;; +infinity, but the following code returns approx 176 + i*pi/4. The
907 ;; reason for the imaginary part is caused by the fact that arg i*y is
908 ;; never 0 since we have positive and negative zeroes.
909
910 (defun complex-atanh (z)
911 "Compute atanh z = (log(1+z) - log(1-z))/2"
912 (declare (number z))
913 (let* (;; Constants
914 (theta #.(/ (sqrt most-positive-double-float) 4.0d0))
915 (rho #.(/ 4.0d0 (sqrt most-positive-double-float)))
916 (half-pi #.(/ pi 2.0d0))
917 (rp (float (realpart z) 1.0d0))
918 (beta (float-sign rp 1.0d0))
919 (x (* beta rp))
920 (y (* beta (- (float (imagpart z) 1.0d0))))
921 (eta 0.0d0)
922 (nu 0.0d0))
923 (declare (double-float theta rho half-pi rp beta y eta nu)
924 (type (double-float 0d0) x))
925 (cond ((or (> x theta)
926 (> (abs y) theta))
927 ;; To avoid overflow...
928 (setf eta (float-sign y half-pi))
929 ;; nu is real part of 1/(x + iy). This is x/(x^2+y^2),
930 ;; which can cause overflow. Arrange this computation so
931 ;; that it won't overflow.
932 (setf nu (let* ((x-bigger (> x (abs y)))
933 (r (if x-bigger (/ y x) (/ x y)))
934 (d (+ 1.0d0 (* r r))))
935 (declare (double-float r d))
936 (if x-bigger
937 (/ (/ x) d)
938 (/ (/ r y) d)))))
939 ((= x 1.0d0)
940 ;; Should this be changed so that if y is zero, eta is set
941 ;; to +infinity instead of approx 176? In any case
942 ;; tanh(176) is 1.0d0 within working precision.
943 (let ((t1 (+ 4d0 (square y)))
944 (t2 (+ (abs y) rho)))
945 (declare (type (double-float 0d0) t1 t2))
946 #+nil
947 (setf eta (log (/ (sqrt (sqrt t1)))
948 (sqrt t2)))
949 (setf eta (* 0.5d0 (log (the (double-float 0.0d0)
950 (/ (sqrt t1) t2)))))
951 (setf nu (* 0.5d0
952 (float-sign y
953 (+ half-pi (atan (* 0.5d0 t2))))))))
954 (t
955 (let ((t1 (+ (abs y) rho)))
956 (declare (double-float t1))
957 ;; Normal case using log1p(x) = log(1 + x)
958 (setf eta (* 0.25d0
959 (%log1p (/ (* 4.0d0 x)
960 (+ (square (- 1.0d0 x))
961 (square t1))))))
962 (setf nu (* 0.5d0
963 (atan (* 2.0d0 y)
964 (- (* (- 1.0d0 x)
965 (+ 1.0d0 x))
966 (square t1))))))))
967 (coerce-to-complex-type (* beta eta)
968 (- (* beta nu))
969 z)))
970
971 (defun complex-tanh (z)
972 "Compute tanh z = sinh z / cosh z"
973 (declare (number z))
974 (let ((x (float (realpart z) 1.0d0))
975 (y (float (imagpart z) 1.0d0)))
976 (declare (double-float x y))
977 (cond ((> (abs x)
978 #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
979 ;; This is more accurate under linux.
980 #+(or linux hpux) #.(/ (+ (%log 2.0d0)
981 (%log most-positive-double-float)) 4d0))
982 (complex (float-sign x)
983 (float-sign y 0.0d0)))
984 (t
985 (let* ((tv (%tan y))
986 (beta (+ 1.0d0 (* tv tv)))
987 (s (sinh x))
988 (rho (sqrt (+ 1.0d0 (* s s)))))
989 (declare (double-float tv s)
990 (type (double-float 0.0d0) beta rho))
991 (if (float-infinity-p (abs tv))
992 (coerce-to-complex-type (/ rho s)
993 (/ tv)
994 z)
995 (let ((den (+ 1.0d0 (* beta s s))))
996 (coerce-to-complex-type (/ (* beta rho s)
997 den)
998 (/ tv den)
999 z))))))))
1000
1001 ;; Kahan says we should only compute the parts needed. Thus, the
1002 ;; realpart's below should only compute the real part, not the whole
1003 ;; complex expression. Doing this can be important because we may get
1004 ;; spurious signals that occur in the part that we are not using.
1005 ;;
1006 ;; However, we take a pragmatic approach and just use the whole
1007 ;; expression.
1008
1009 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1010 ;; it's the conjugate of the square root or the square root of the
1011 ;; conjugate. This needs to be checked.
1012
1013 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1014 ;; same as (sqrt (conjugate z)) for all z. This follows because
1015 ;;
1016 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1017 ;;
1018 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1019 ;;
1020 ;;.and these two expressions are equal if and only if arg conj z =
1021 ;;-arg z, which is clearly true for all z.
1022
1023 (defun complex-acos (z)
1024 "Compute acos z = pi/2 - asin z
1025
1026 Z may be any number, but the result is always a complex."
1027 (declare (number z))
1028 (let ((sqrt-1+z (complex-sqrt (+ 1 z)))
1029 (sqrt-1-z (complex-sqrt (- 1 z))))
1030 (with-float-traps-masked (:divide-by-zero)
1031 (complex (* 2 (atan (/ (realpart sqrt-1-z)
1032 (realpart sqrt-1+z))))
1033 (asinh (imagpart (* (conjugate sqrt-1+z)
1034 sqrt-1-z)))))))
1035
1036 (defun complex-acosh (z)
1037 "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1038
1039 Z may be any number, but the result is always a complex."
1040 (declare (number z))
1041 (let ((sqrt-z-1 (complex-sqrt (- z 1)))
1042 (sqrt-z+1 (complex-sqrt (+ z 1))))
1043 (with-float-traps-masked (:divide-by-zero)
1044 (complex (asinh (realpart (* (conjugate sqrt-z-1)
1045 sqrt-z+1)))
1046 (* 2 (atan (/ (imagpart sqrt-z-1)
1047 (realpart sqrt-z+1))))))))
1048
1049
1050 (defun complex-asin (z)
1051 "Compute asin z = asinh(i*z)/i
1052
1053 Z may be any number, but the result is always a complex."
1054 (declare (number z))
1055 (let ((sqrt-1-z (complex-sqrt (- 1 z)))
1056 (sqrt-1+z (complex-sqrt (+ 1 z))))
1057 (with-float-traps-masked (:divide-by-zero)
1058 (complex (atan (/ (realpart z)
1059 (realpart (* sqrt-1-z sqrt-1+z))))
1060 (asinh (imagpart (* (conjugate sqrt-1-z)
1061 sqrt-1+z)))))))
1062
1063 (defun complex-asinh (z)
1064 "Compute asinh z = log(z + sqrt(1 + z*z))
1065
1066 Z may be any number, but the result is always a complex."
1067 (declare (number z))
1068 ;; asinh z = -i * asin (i*z)
1069 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1070 (result (complex-asin iz)))
1071 (complex (imagpart result)
1072 (- (realpart result)))))
1073
1074 (defun complex-atan (z)
1075 "Compute atan z = atanh (i*z) / i
1076
1077 Z may be any number, but the result is always a complex."
1078 (declare (number z))
1079 ;; atan z = -i * atanh (i*z)
1080 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1081 (result (complex-atanh iz)))
1082 (complex (imagpart result)
1083 (- (realpart result)))))
1084
1085 (defun complex-tan (z)
1086 "Compute tan z = -i * tanh(i * z)
1087
1088 Z may be any number, but the result is always a complex."
1089 (declare (number z))
1090 ;; tan z = -i * tanh(i*z)
1091 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1092 (result (complex-tanh iz)))
1093 (complex (imagpart result)
1094 (- (realpart result)))))
1095 )

  ViewVC Help
Powered by ViewVC 1.1.5