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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.19.2.2 - (show annotations)
Mon Sep 8 00:22:00 1997 UTC (16 years, 7 months ago) by dtc
Branch: RELENG_18
CVS Tags: RELEASE_18a
Changes since 1.19.2.1: +3 -3 lines
Check the accrued exceptions (sticky-bits) within
with-float-traps-masked; more reliable than the exceptions-byte on
many ports.
Call %tan directly in complex-tanh, saves a compiler note on x86.
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.19.2.2 1997/09/08 00:22:00 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 ((real-expt (base power rtype)
205 (let* ((fbase (coerce base 'double-float))
206 (fpower (coerce power 'double-float))
207 (res (coerce (%pow fbase fpower) rtype)))
208 (if (and (zerop res) (minusp fbase))
209 (multiple-value-bind (re im)
210 (complex-pow fbase fpower)
211 (%make-complex (coerce re rtype) (coerce im rtype)))
212 res)))
213 (complex-pow (fbase fpower)
214 (let ((pow (%pow (- fbase) fpower))
215 (fpower*pi (* fpower pi)))
216 (values (* pow (%cos fpower*pi))
217 (* pow (%sin fpower*pi))))))
218 (declare (inline real-expt))
219 (number-dispatch ((base number) (power number))
220 (((foreach fixnum (or bignum ratio) (complex rational)) integer)
221 (intexp base power))
222 (((foreach single-float double-float) rational)
223 (real-expt base power '(dispatch-type base)))
224 (((foreach fixnum (or bignum ratio) single-float)
225 (foreach ratio single-float))
226 (real-expt base power 'single-float))
227 (((foreach fixnum (or bignum ratio) single-float double-float)
228 double-float)
229 (real-expt base power 'double-float))
230 ((double-float single-float)
231 (real-expt base power 'double-float))
232 (((foreach (complex rational) (complex float)) rational)
233 (* (expt (abs base) power)
234 (cis (* power (phase base)))))
235 (((foreach fixnum (or bignum ratio) single-float double-float)
236 complex)
237 (exp (* power (log base))))
238 (((foreach (complex float) (complex rational)) number)
239 (exp (* power (log base))))))))
240
241 (defun log (number &optional (base nil base-p))
242 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
243 (if base-p
244 (if (zerop base)
245 base ; ANSI spec
246 (/ (log number) (log base)))
247 (number-dispatch ((number number))
248 (((foreach fixnum bignum ratio))
249 (if (minusp number)
250 (complex (log (- number)) (coerce pi 'single-float))
251 (coerce (%log (coerce number 'double-float)) 'single-float)))
252 (((foreach single-float double-float))
253 ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
254 ;; Since this doesn't seem to be an implementation issue
255 ;; I (pw) take the Kahan result.
256 (if (< (float-sign number)
257 (coerce 0 '(dispatch-type number)))
258 (complex (log (- number)) (coerce pi '(dispatch-type number)))
259 (coerce (%log (coerce number 'double-float))
260 '(dispatch-type number))))
261 ((complex)
262 (complex-log number)))))
263
264 (defun sqrt (number)
265 "Return the square root of NUMBER."
266 (number-dispatch ((number number))
267 (((foreach fixnum bignum ratio))
268 (if (minusp number)
269 (complex-sqrt number)
270 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
271 (((foreach single-float double-float))
272 ;; NOTE there is a problem with (at least x86 NPX) of what result
273 ;; should be returned for (sqrt -0.0). The x86 hardware FSQRT
274 ;; instruction returns -0d0. The result is that Python will perhaps
275 ;; note the following test in generic sqrt, non-negatively constrained
276 ;; float types will be passed to FSQRT (or libm on other boxes).
277 ;; So, in the interest of consistency of compiled and interpreted
278 ;; codes, the following test is disabled for now. Maybe the float-sign
279 ;; test could be moved to the optimization codes.
280 (if (< (#+nil float-sign #-nil identity number)
281 (coerce 0 '(dispatch-type number)))
282 (complex-sqrt number)
283 (coerce (%sqrt (coerce number 'double-float))
284 '(dispatch-type number))))
285 ((complex)
286 (complex-sqrt number))))
287
288
289 ;;;; Trigonometic and Related Functions
290
291 (defun abs (number)
292 "Returns the absolute value of the number."
293 (number-dispatch ((number number))
294 (((foreach single-float double-float fixnum rational))
295 (abs number))
296 ((complex)
297 (let ((rx (realpart number))
298 (ix (imagpart number)))
299 (etypecase rx
300 (rational
301 (sqrt (+ (* rx rx) (* ix ix))))
302 (single-float
303 (coerce (%hypot (coerce rx 'double-float)
304 (coerce ix 'double-float))
305 'single-float))
306 (double-float
307 (%hypot rx ix)))))))
308
309 (defun phase (number)
310 "Returns the angle part of the polar representation of a complex number.
311 For complex numbers, this is (atan (imagpart number) (realpart number)).
312 For non-complex positive numbers, this is 0. For non-complex negative
313 numbers this is PI."
314 (etypecase number
315 (rational
316 (if (minusp number)
317 (coerce pi 'single-float)
318 0.0f0))
319 (single-float
320 (if (minusp (float-sign number))
321 (coerce pi 'single-float)
322 0.0f0))
323 (double-float
324 (if (minusp (float-sign number))
325 (coerce pi 'double-float)
326 0.0d0))
327 (complex
328 (atan (imagpart number) (realpart number)))))
329
330
331 (defun sin (number)
332 "Return the sine of NUMBER."
333 (number-dispatch ((number number))
334 (handle-reals %sin number)
335 ((complex)
336 (let ((x (realpart number))
337 (y (imagpart number)))
338 (complex (* (sin x) (cosh y))
339 (* (cos x) (sinh y)))))))
340
341 (defun cos (number)
342 "Return the cosine of NUMBER."
343 (number-dispatch ((number number))
344 (handle-reals %cos number)
345 ((complex)
346 (let ((x (realpart number))
347 (y (imagpart number)))
348 (complex (* (cos x) (cosh y))
349 (- (* (sin x) (sinh y))))))))
350
351 (defun tan (number)
352 "Return the tangent of NUMBER."
353 (number-dispatch ((number number))
354 (handle-reals %tan number)
355 ((complex)
356 (complex-tan number))))
357
358 (defun cis (theta)
359 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
360 (if (complexp theta)
361 (error "Argument to CIS is complex: ~S" theta)
362 (complex (cos theta) (sin theta))))
363
364 (defun asin (number)
365 "Return the arc sine of NUMBER."
366 (number-dispatch ((number number))
367 ((rational)
368 (if (or (> number 1) (< number -1))
369 (complex-asin number)
370 (coerce (%asin (coerce number 'double-float)) 'single-float)))
371 (((foreach single-float double-float))
372 (if (or (> number (coerce 1 '(dispatch-type number)))
373 (< number (coerce -1 '(dispatch-type number))))
374 (complex-asin number)
375 (coerce (%asin (coerce number 'double-float))
376 '(dispatch-type number))))
377 ((complex)
378 (complex-asin number))))
379
380 (defun acos (number)
381 "Return the arc cosine of NUMBER."
382 (number-dispatch ((number number))
383 ((rational)
384 (if (or (> number 1) (< number -1))
385 (complex-acos number)
386 (coerce (%acos (coerce number 'double-float)) 'single-float)))
387 (((foreach single-float double-float))
388 (if (or (> number (coerce 1 '(dispatch-type number)))
389 (< number (coerce -1 '(dispatch-type number))))
390 (complex-acos number)
391 (coerce (%acos (coerce number 'double-float))
392 '(dispatch-type number))))
393 ((complex)
394 (complex-acos number))))
395
396
397 (defun atan (y &optional (x nil xp))
398 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
399 (if xp
400 (flet ((atan2 (y x)
401 (declare (type double-float y x)
402 (values double-float))
403 (if (zerop x)
404 (if (zerop y)
405 (if (plusp (float-sign x))
406 y
407 (float-sign y pi))
408 (float-sign y (/ pi 2)))
409 (%atan2 y x))))
410 (number-dispatch ((y number) (x number))
411 ((double-float
412 (foreach double-float single-float fixnum bignum ratio))
413 (atan2 y (coerce x 'double-float)))
414 (((foreach single-float fixnum bignum ratio)
415 double-float)
416 (atan2 (coerce y 'double-float) x))
417 (((foreach single-float fixnum bignum ratio)
418 (foreach single-float fixnum bignum ratio))
419 (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
420 'single-float))))
421 (number-dispatch ((y number))
422 (handle-reals %atan y)
423 ((complex)
424 (complex-atan y)))))
425
426 ;; It seems that everyone has a C version of sinh, cosh, and
427 ;; tanh. Let's use these for reals because the original
428 ;; implementations based on the definitions lose big in round-off
429 ;; error. These bad definitions also mean that sin and cos for
430 ;; complex numbers can also lose big.
431
432 #+nil
433 (defun sinh (number)
434 "Return the hyperbolic sine of NUMBER."
435 (/ (- (exp number) (exp (- number))) 2))
436
437 (defun sinh (number)
438 "Return the hyperbolic sine of NUMBER."
439 (number-dispatch ((number number))
440 (handle-reals %sinh number)
441 ((complex)
442 (let ((x (realpart number))
443 (y (imagpart number)))
444 (complex (* (sinh x) (cos y))
445 (* (cosh x) (sin y)))))))
446
447 #+nil
448 (defun cosh (number)
449 "Return the hyperbolic cosine of NUMBER."
450 (/ (+ (exp number) (exp (- number))) 2))
451
452 (defun cosh (number)
453 "Return the hyperbolic cosine of NUMBER."
454 (number-dispatch ((number number))
455 (handle-reals %cosh number)
456 ((complex)
457 (let ((x (realpart number))
458 (y (imagpart number)))
459 (complex (* (cosh x) (cos y))
460 (* (sinh x) (sin y)))))))
461
462 (defun tanh (number)
463 "Return the hyperbolic tangent of NUMBER."
464 (number-dispatch ((number number))
465 (handle-reals %tanh number)
466 ((complex)
467 (complex-tanh number))))
468
469 (defun asinh (number)
470 "Return the hyperbolic arc sine of NUMBER."
471 (number-dispatch ((number number))
472 (handle-reals %asinh number)
473 ((complex)
474 (complex-asinh number))))
475
476 (defun acosh (number)
477 "Return the hyperbolic arc cosine of NUMBER."
478 (number-dispatch ((number number))
479 ((rational)
480 ;; acosh is complex if number < 1
481 (if (< number 1)
482 (complex-acosh number)
483 (coerce (%acosh (coerce number 'double-float)) 'single-float)))
484 (((foreach single-float double-float))
485 (if (< number (coerce 1 '(dispatch-type number)))
486 (complex-acosh number)
487 (coerce (%acosh (coerce number 'double-float))
488 '(dispatch-type number))))
489 ((complex)
490 (complex-acosh number))))
491
492 (defun atanh (number)
493 "Return the hyperbolic arc tangent of NUMBER."
494 (number-dispatch ((number number))
495 ((rational)
496 ;; atanh is complex if |number| > 1
497 (if (or (> number 1) (< number -1))
498 (complex-atanh number)
499 (coerce (%atanh (coerce number 'double-float)) 'single-float)))
500 (((foreach single-float double-float))
501 (if (or (> number (coerce 1 '(dispatch-type number)))
502 (< number (coerce -1 '(dispatch-type number))))
503 (complex-atanh number)
504 (coerce (%atanh (coerce number 'double-float))
505 '(dispatch-type number))))
506 ((complex)
507 (complex-atanh number))))
508
509 ;;; HP-UX does not supply a C version of log1p, so
510 ;;; use the definition.
511
512 #+hpux
513 (declaim (inline %log1p))
514 #+hpux
515 (defun %log1p (number)
516 (declare (double-float number)
517 (optimize (speed 3) (safety 0)))
518 (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
519
520
521 #+old-elfun
522 (progn
523 ;;; Here are the old definitions of the special functions, for
524 ;;; complex-valued arguments. Some of these functions suffer from
525 ;;; severe round-off error or unnecessary overflow.
526
527 (proclaim '(inline mult-by-i))
528 (defun mult-by-i (number)
529 (complex (- (imagpart number))
530 (realpart number)))
531
532 (defun complex-sqrt (x)
533 (exp (/ (log x) 2)))
534
535 (defun complex-log (x)
536 (complex (log (abs x))
537 (phase x)))
538
539 (defun complex-atanh (number)
540 (/ (- (log (1+ number)) (log (- 1 number))) 2))
541
542 (defun complex-tanh (number)
543 (/ (- (exp number) (exp (- number)))
544 (+ (exp number) (exp (- number)))))
545
546 (defun complex-acos (number)
547 (* -2 (mult-by-i (log (+ (sqrt (/ (1+ number) 2))
548 (mult-by-i (sqrt (/ (- 1 number) 2))))))))
549
550 (defun complex-acosh (number)
551 (* 2 (log (+ (sqrt (/ (1+ number) 2)) (sqrt (/ (1- number) 2))))))
552
553 (defun complex-asin (number)
554 (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
555
556 (defun complex-asinh (number)
557 (log (+ number (sqrt (1+ (* number number))))))
558
559 (defun complex-atan (y)
560 (let ((im (imagpart y))
561 (re (realpart y)))
562 (/ (- (log (complex (- 1 im) re))
563 (log (complex (+ 1 im) (- re))))
564 (complex 0 2))))
565
566 (defun complex-tan (number)
567 (let* ((num (sin number))
568 (denom (cos number)))
569 (if (zerop denom) (error "~S undefined tangent." number)
570 (/ num denom))))
571 )
572
573 #-old-specfun
574 (progn
575 ;;;;
576 ;;;; This is a set of routines that implement many elementary
577 ;;;; transcendental functions as specified by ANSI Common Lisp. The
578 ;;;; implementation is based on Kahan's paper.
579 ;;;;
580 ;;;; I believe I have accurately implemented the routines and are
581 ;;;; correct, but you may want to check for your self.
582 ;;;;
583 ;;;; These functions are written for CMU Lisp and take advantage of
584 ;;;; some of the features available there. It may be possible,
585 ;;;; however, to port this to other Lisps.
586 ;;;;
587 ;;;; Some functions are significantly more accurate than the original
588 ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
589 ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
590 ;;;; answer is pi + i*log(2-sqrt(3)).
591 ;;;;
592 ;;;; All of the implemented functions will take any number for an
593 ;;;; input, but the result will always be a either a complex
594 ;;;; single-float or a complex double-float.
595 ;;;;
596 ;;;; General functions
597 ;;;; complex-sqrt
598 ;;;; complex-log
599 ;;;; complex-atanh
600 ;;;; complex-tanh
601 ;;;; complex-acos
602 ;;;; complex-acosh
603 ;;;; complex-asin
604 ;;;; complex-asinh
605 ;;;; complex-atan
606 ;;;; complex-tan
607 ;;;;
608 ;;;; Utility functions:
609 ;;;; scalb logb
610 ;;;;
611 ;;;; Internal functions:
612 ;;;; square coerce-to-complex-type cssqs complex-log-scaled
613 ;;;;
614 ;;;;
615 ;;;; Please send any bug reports, comments, or improvements to Raymond
616 ;;;; Toy at toy@rtp.ericsson.se.
617 ;;;;
618 ;;;; References
619 ;;;;
620 ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
621 ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
622 ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
623 ;;;; Press, 1987
624 ;;;;
625 (declaim (inline square))
626 (declaim (ftype (function (double-float) (double-float 0d0)) square))
627 (defun square (x)
628 (declare (double-float x)
629 (values (double-float 0d0)))
630 (* x x))
631
632 ;; If you have these functions in libm, perhaps they should be used
633 ;; instead of these Lisp versions. These versions are probably good
634 ;; enough, especially since they are portable.
635
636 (declaim (inline scalb))
637 (defun scalb (x n)
638 "Compute 2^N * X without compute 2^N first (use properties of the
639 underlying floating-point format"
640 (declare (type double-float x)
641 (type double-float-exponent n))
642 (scale-float x n))
643
644 (defun logb (x)
645 "Compute an integer N such that 1 <= |2^N * x| < 2.
646 For the special cases, the following values are used:
647
648 x logb
649 NaN NaN
650 +/- infinity +infinity
651 0 -infinity
652 "
653 (declare (type double-float x))
654 (cond ((float-trapping-nan-p x)
655 x)
656 ((float-infinity-p x)
657 #.ext:double-float-positive-infinity)
658 ((zerop x)
659 ;; The answer is negative infinity, but we are supposed to
660 ;; signal divide-by-zero.
661 ;; (error 'division-by-zero :operation 'logb :operands (list x))
662 (/ -1.0d0 x)
663 )
664 (t
665 (multiple-value-bind (signif expon sign)
666 (decode-float x)
667 (declare (ignore signif sign))
668 ;; decode-float is almost right, except that the exponent
669 ;; is off by one
670 (1- expon)))))
671
672 ;; This function is used to create a complex number of the appropriate
673 ;; type.
674
675 (declaim (inline coerce-to-complex-type))
676 (defun coerce-to-complex-type (x y z)
677 "Create complex number with real part X and imaginary part Y such that
678 it has the same type as Z. If Z has type (complex rational), the X
679 and Y are coerced to single-float."
680 (declare (double-float x y)
681 (number z))
682 (if (subtypep (type-of (realpart z)) 'double-float)
683 (complex x y)
684 ;; Convert anything that's not a double-float to a single-float.
685 (complex (float x 1.0)
686 (float y 1.0))))
687
688 (defun cssqs (z)
689 ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
690 ;; result is r + i*k, where k is an integer.
691
692 ;; Save all FP flags
693 (let ((x (float (realpart z) 1d0))
694 (y (float (imagpart z) 1d0))
695 (k 0)
696 (rho 0d0))
697 (declare (double-float x y)
698 (type (double-float 0d0) rho)
699 (fixnum k))
700 ;; Would this be better handled using an exception handler to
701 ;; catch the overflow or underflow signal? For now, we turn all
702 ;; traps off and look at the accrued exceptions to see if any
703 ;; signal would have been raised.
704 (with-float-traps-masked (:underflow :overflow)
705 (setf rho (+ (square x) (square y)))
706 (cond ((and (or (float-trapping-nan-p rho)
707 (float-infinity-p rho))
708 (or (float-infinity-p (abs x))
709 (float-infinity-p (abs y))))
710 (setf rho #.ext:double-float-positive-infinity))
711 ((let ((threshold #.(/ least-positive-double-float
712 double-float-epsilon))
713 (traps (ldb vm::float-sticky-bits
714 (vm:floating-point-modes))))
715 ;; Overflow raised or (underflow raised and rho <
716 ;; lambda/eps)
717 (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
718 (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
719 (< rho threshold))))
720 (setf k (logb (max (abs x) (abs y))))
721 (setf rho (+ (square (scalb x (- k)))
722 (square (scalb y (- k))))))))
723 (values rho k)))
724
725 (defun complex-sqrt (z)
726 "Principle square root of Z
727
728 Z may be any number, but the result is always a complex."
729 (declare (number z))
730 (multiple-value-bind (rho k)
731 (cssqs z)
732 (declare (type (double-float 0d0) rho)
733 (fixnum k))
734 (let ((x (float (realpart z) 1.0d0))
735 (y (float (imagpart z) 1.0d0))
736 (eta 0d0)
737 (nu 0d0))
738 (declare (double-float x y eta nu))
739
740 (if (not (float-trapping-nan-p x))
741 (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
742
743 (cond ((oddp k)
744 (setf k (ash k -1)))
745 (t
746 (setf k (1- (ash k -1)))
747 (setf rho (+ rho rho))))
748
749 (setf rho (scalb (sqrt rho) k))
750
751 (setf eta rho)
752 (setf nu y)
753
754 (when (/= rho 0d0)
755 (when (not (float-infinity-p (abs nu)))
756 (setf nu (/ (/ nu rho) 2d0)))
757 (when (< x 0d0)
758 (setf eta (abs nu))
759 (setf nu (float-sign y rho))))
760 (coerce-to-complex-type eta nu z))))
761
762 (defun complex-log-scaled (z j)
763 "Compute log(2^j*z).
764
765 This is for use with J /= 0 only when |z| is huge."
766 (declare (number z)
767 (fixnum j))
768 ;; The constants t0, t1, t2 should be evaluated to machine
769 ;; precision. In addition, Kahan says the accuracy of log1p
770 ;; influences the choices of these constants but doesn't say how to
771 ;; choose them. We'll just assume his choices matches our
772 ;; implementation of log1p.
773 (let ((t0 #.(/ 1 (sqrt 2.0d0)))
774 (t1 1.2d0)
775 (t2 3d0)
776 (ln2 #.(log 2d0))
777 (x (float (realpart z) 1.0d0))
778 (y (float (imagpart z) 1.0d0)))
779 (multiple-value-bind (rho k)
780 (cssqs z)
781 (declare (type (double-float 0d0) rho)
782 (fixnum k))
783 (let ((beta (max (abs x) (abs y)))
784 (theta (min (abs x) (abs y))))
785 (declare (type (double-float 0d0) beta theta))
786 (if (and (zerop k)
787 (< t0 beta)
788 (or (<= beta t1)
789 (< rho t2)))
790 (setf rho (/ (%log1p (+ (* (- beta 1.0d0)
791 (+ beta 1.0d0))
792 (* theta theta)))
793 2d0))
794 (setf rho (+ (/ (log rho) 2d0)
795 (* (+ k j) ln2))))
796 (setf theta (atan y x))
797 (coerce-to-complex-type rho theta z)))))
798
799 (defun complex-log (z)
800 "Log of Z = log |Z| + i * arg Z
801
802 Z may be any number, but the result is always a complex."
803 (declare (number z))
804 (complex-log-scaled z 0))
805
806 ;; Let us note the following "strange" behavior. atanh 1.0d0 is
807 ;; +infinity, but the following code returns approx 176 + i*pi/4. The
808 ;; reason for the imaginary part is caused by the fact that arg i*y is
809 ;; never 0 since we have positive and negative zeroes.
810
811 (defun complex-atanh (z)
812 "Compute atanh z = (log(1+z) - log(1-z))/2"
813 (declare (number z))
814 (let* (;; Constants
815 (theta #.(/ (sqrt most-positive-double-float) 4.0d0))
816 (rho #.(/ 4.0d0 (sqrt most-positive-double-float)))
817 (half-pi #.(/ pi 2.0d0))
818 (rp (float (realpart z) 1.0d0))
819 (beta (float-sign rp 1.0d0))
820 (x (* beta rp))
821 (y (* beta (- (float (imagpart z) 1.0d0))))
822 (eta 0.0d0)
823 (nu 0.0d0))
824 (declare (double-float theta rho half-pi rp beta y eta nu)
825 (type (double-float 0d0) x))
826 (cond ((or (> x theta)
827 (> (abs y) theta))
828 ;; To avoid overflow...
829 (setf eta (float-sign y half-pi))
830 ;; nu is real part of 1/(x + iy). This is x/(x^2+y^2),
831 ;; which can cause overflow. Arrange this computation so
832 ;; that it won't overflow.
833 (setf nu (let* ((x-bigger (> x (abs y)))
834 (r (if x-bigger (/ y x) (/ x y)))
835 (d (+ 1.0d0 (* r r))))
836 (declare (double-float r d))
837 (if x-bigger
838 (/ (/ x) d)
839 (/ (/ r y) d)))))
840 ((= x 1.0d0)
841 ;; Should this be changed so that if y is zero, eta is set
842 ;; to +infinity instead of approx 176? In any case
843 ;; tanh(176) is 1.0d0 within working precision.
844 (let ((t1 (+ 4d0 (square y)))
845 (t2 (+ (abs y) rho)))
846 (declare (type (double-float 0d0) t1 t2))
847 #+nil
848 (setf eta (log (/ (sqrt (sqrt t1)))
849 (sqrt t2)))
850 (setf eta (* 0.5d0 (log (the (double-float 0.0d0)
851 (/ (sqrt t1) t2)))))
852 (setf nu (* 0.5d0
853 (float-sign y
854 (+ half-pi (atan (* 0.5d0 t2))))))))
855 (t
856 (let ((t1 (+ (abs y) rho)))
857 (declare (double-float t1))
858 ;; Normal case using log1p(x) = log(1 + x)
859 (setf eta (* 0.25d0
860 (%log1p (/ (* 4.0d0 x)
861 (+ (square (- 1.0d0 x))
862 (square t1))))))
863 (setf nu (* 0.5d0
864 (atan (* 2.0d0 y)
865 (- (* (- 1.0d0 x)
866 (+ 1.0d0 x))
867 (square t1))))))))
868 (coerce-to-complex-type (* beta eta)
869 (- (* beta nu))
870 z)))
871
872 (defun complex-tanh (z)
873 "Compute tanh z = sinh z / cosh z"
874 (declare (number z))
875 (let ((x (float (realpart z) 1.0d0))
876 (y (float (imagpart z) 1.0d0)))
877 (declare (double-float x y))
878 (cond ((> (abs x)
879 #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
880 ;; This is more accurate under linux.
881 #+(or linux hpux) #.(/ (+ (%log 2.0d0)
882 (%log most-positive-double-float)) 4d0))
883 (complex (float-sign x)
884 (float-sign y 0.0d0)))
885 (t
886 (let* ((tv (%tan y))
887 (beta (+ 1.0d0 (* tv tv)))
888 (s (sinh x))
889 (rho (sqrt (+ 1.0d0 (* s s)))))
890 (declare (double-float tv s)
891 (type (double-float 0.0d0) beta rho))
892 (if (float-infinity-p (abs tv))
893 (coerce-to-complex-type (/ rho s)
894 (/ tv)
895 z)
896 (let ((den (+ 1.0d0 (* beta s s))))
897 (coerce-to-complex-type (/ (* beta rho s)
898 den)
899 (/ tv den)
900 z))))))))
901
902 ;; Kahan says we should only compute the parts needed. Thus, the
903 ;; realpart's below should only compute the real part, not the whole
904 ;; complex expression. Doing this can be important because we may get
905 ;; spurious signals that occur in the part that we are not using.
906 ;;
907 ;; However, we take a pragmatic approach and just use the whole
908 ;; expression.
909
910 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
911 ;; it's the conjugate of the square root or the square root of the
912 ;; conjugate. This needs to be checked.
913
914 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
915 ;; same as (sqrt (conjugate z)) for all z. This follows because
916 ;;
917 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
918 ;;
919 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
920 ;;
921 ;;.and these two expressions are equal if and only if arg conj z =
922 ;;-arg z, which is clearly true for all z.
923
924 (defun complex-acos (z)
925 "Compute acos z = pi/2 - asin z
926
927 Z may be any number, but the result is always a complex."
928 (declare (number z))
929 (let ((sqrt-1+z (complex-sqrt (+ 1 z)))
930 (sqrt-1-z (complex-sqrt (- 1 z))))
931 (with-float-traps-masked (:divide-by-zero)
932 (complex (* 2 (atan (/ (realpart sqrt-1-z)
933 (realpart sqrt-1+z))))
934 (asinh (imagpart (* (conjugate sqrt-1+z)
935 sqrt-1-z)))))))
936
937 (defun complex-acosh (z)
938 "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
939
940 Z may be any number, but the result is always a complex."
941 (declare (number z))
942 (let ((sqrt-z-1 (complex-sqrt (- z 1)))
943 (sqrt-z+1 (complex-sqrt (+ z 1))))
944 (with-float-traps-masked (:divide-by-zero)
945 (complex (asinh (realpart (* (conjugate sqrt-z-1)
946 sqrt-z+1)))
947 (* 2 (atan (/ (imagpart sqrt-z-1)
948 (realpart sqrt-z+1))))))))
949
950
951 (defun complex-asin (z)
952 "Compute asin z = asinh(i*z)/i
953
954 Z may be any number, but the result is always a complex."
955 (declare (number z))
956 (let ((sqrt-1-z (complex-sqrt (- 1 z)))
957 (sqrt-1+z (complex-sqrt (+ 1 z))))
958 (with-float-traps-masked (:divide-by-zero)
959 (complex (atan (/ (realpart z)
960 (realpart (* sqrt-1-z sqrt-1+z))))
961 (asinh (imagpart (* (conjugate sqrt-1-z)
962 sqrt-1+z)))))))
963
964 (defun complex-asinh (z)
965 "Compute asinh z = log(z + sqrt(1 + z*z))
966
967 Z may be any number, but the result is always a complex."
968 (declare (number z))
969 ;; asinh z = -i * asin (i*z)
970 (let* ((iz (complex (- (imagpart z)) (realpart z)))
971 (result (complex-asin iz)))
972 (complex (imagpart result)
973 (- (realpart result)))))
974
975 (defun complex-atan (z)
976 "Compute atan z = atanh (i*z) / i
977
978 Z may be any number, but the result is always a complex."
979 (declare (number z))
980 ;; atan z = -i * atanh (i*z)
981 (let* ((iz (complex (- (imagpart z)) (realpart z)))
982 (result (complex-atanh iz)))
983 (complex (imagpart result)
984 (- (realpart result)))))
985
986 (defun complex-tan (z)
987 "Compute tan z = -i * tanh(i * z)
988
989 Z may be any number, but the result is always a complex."
990 (declare (number z))
991 ;; tan z = -i * tanh(i*z)
992 (let* ((iz (complex (- (imagpart z)) (realpart z)))
993 (result (complex-tanh iz)))
994 (complex (imagpart result)
995 (- (realpart result)))))
996 )

  ViewVC Help
Powered by ViewVC 1.1.5